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SUWARY 


This investigation combines an incompressible inviscid and a viscous flow 
calculation procedure by assuming the viscous correction to the inviscid flow 
pressure distribution is small (weak interaction) to predict the flow about 
airfoils oscillating in pitch and heave. The calculations made in the 
investigation predict the detailed viscous flow regions including transition 
and separation phenomena and provide a detailed analysis of leading edge 
separation, transition, and reattachment. Results from the calculation show 
the leading edge viscous flew field to be quasi -steady altnough the imposed 
inviscid pressure distribution shows significant unsteady effects. Although 
unable to predict the flow field about a stalled airfoil, the indications 
are that the present procedure can indicate the onset of catastrophic flow 
separation. 


INTRODUCTION 


The phenomenon of dynamic stall, which is largely controlled by the 
viscous boundary layer in direct contact with the airfoil surface, plays an 
important role in the successful design and operation of helicopter rotor 
blades. Under high-speed flight conditions the retreating rotor blades are 
subject to a diminished dynamic pressure and, as a result, high blade 
performance requires large lift coefficients to be present in the retreating 
port i cxi of the rotor disc. Ihese large lift coefficients are generated 
through large incidence angles often exceeding the maximum angle for which the 
boundary layer can remain completely attached to the airfoil surface even under 
dynamic conditions. When a significant amount of boundary layer separation 
appears, the airfoil experiences a deterioration in performance which is termed 
stall. Stall is most easily described in terms of a lift coefficient -incidence 


♦The contract research effort which has lead to the results in this report 
was financially supported by USAAMRDL (Langley Directorate) 



angle relation. At low incidence angles no significant amount of boundary 
layer separation is present and the lift coefficient varies linearly with 
incidence, as is predicted by inviscid flow theory. At some given incidence 
the lift coefficient -incidence curve becomes nonlinear as the lift increases 
less quickly than incidence; this decrease in the lift -incidence slope is 
accompanied by a thickening of the viscous boundary layer and perhaps even by 
boundary layer separation which causes the lift to vary from its inviscid 
value. Further increases in incidence lead to larger decrements in the lift 
coefficient from the linear lift -incidence relation of potential flow and 
eventually an incidence corresponding to a maximum lift coefficient is 
reached; after this maximum is reached any further increase in incidence is 
accompanied by a decrease in lift. At these higher incidence angles the 
viscous flew about the airfoil is characterized by large separated regions 
along the suction side of the airfoil and in the airfoil wake, clearly 
indicating a relationship between boundary layer separation and airfoil stall. 
In addition to the behavior of the lift coefficient during stall, the moment 
coefficient about the quarter chord point shows a large change from its nearly 
zero value characteristic of unstalled flow, indicating a significant shift in 
the center of pressure. 

The performance of the airfoil during dynamic stall plays an important 
role in determining the overall helicopter performance. Obviously, the lift 
is highly-dependent upon airfoil peiformance during stall and, furthermore, 
blade fatigue stress, blade flutter, and aircraft vibration are significantly 
affected by the periodic aerodynamic loading and unloading as the blade 
proceeds about the rotor disc. Thus, an accurate procedure for predicting 
the unsteady flew about an airfoil during dynamic stall would represent a 
significant input to a rotor design system. 

McCullough and Gault (ref. l) have postulated three types of stall for 
airfoils in steady flow; these are leading-edge stall, trailing-edge stall, 
and thin-airfoil stall. The first of these, leading -edge stall, is supposedly 
related to the formation of a separation bubble in the vicinity of the airfoil 
leading edge. For leading-edge stall it is conjectured by McCullough and 
Gault (ref. l) that as incidence increases the bubble moves upstream until an 
incidence is reached at which the bubble suddenly bursts and the flow separates 
from the airfoil surface. The bursting process is accompanied by a sudden 
loss in lift and decrease in airfoil performance. In contrast to leading -edge 
stall, which is supposedly associated with separation at the leading edge of 
the airfoil, trailing-edge stall, which usually occurs on relatively thick 
airfoils, is associated with the separation of the boundary layer on the aft 
portion of the airfoil. Under most operating conditions, trailing-edge stall 
is associated with the separation of a turbulent rather than a laminar 
boundary layer. At low incidence no trailing-edge separation occurs. However, 
at some given incidence the boundary layer separates in the vicinity of the 
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trailing edge and, if the incidence continues to increase, the separation 
point moves forward causing a gradual decrement in airfoil performance. The 
effects associated with trailing -edge stall are considerably more gradual than 
those associated with leading-edge stall. The final type of stall conjectured 
by McCullough and Gault (ref. 1) is thin-airfoil stall, which, as the name 
implies, occurs on thin airfoils. Like leading-edge stall, thin-airfoil stall 
is associated with a separation bubble in the vicinity of the leading edge. 
However, ref. 1 suggests that in the case of thin-airfoil stall the bubble 
grows with increasing incidence angle, whereas for leading-edge stall the 
bubble moves upstream and can even shorten with increasing incidence angle. 

The three different types of stall discussed in ref. 1 appear to be 
associated with three different types of viscous separation . However, boundary 
layer separation is an extremely sensitive phenomenon, the nature of which can 
be significantly altered by changes in the applied pressure distribution, free- 
stream turbulence level, Reynolds number, etc., and, therefore, it is not 
unreasonable to expect an airfoil which exhibits om type of stall under a 
given set of conditions to exhibit a different type of stall under different 
conditions. It may even be possible for an airfoil to be subject to leading- 
edge and trailing-edge stall simultaneously. Obviously, since stall is 
heavily dependent upon the extremely sensitive viscous separation mechanism, 
it is questionable how well mechanisms of stall deduced from any specific set 
of data correspond to the mechanisms of stall under different operating 
conditions . 

Although the suggestions of McCullough and Gault (ref. l) apply to static 
stall, it seems reasonable to suppose* that some of these same mechanisms are 
present in dynamic stall. However, important differences do exist between 
the static and dynamic cases (e.g., refs. 2 and 3 ). First of all, the maximum 
incidence angle which the airfoil can tolerate before the linear lift 
coefficient -incidence angle relation breaks down is significantly higher in 
the dynamic case than in the static case, indicating that a delay in stall is 
obtained through dynamic phenomena. The maximum lift coefficient obtained 
under dynamic conditions in general is significantly greater than that obtained 
under static conditions. In addition, under dynamic conditions, stall shews 
a definite hystersis effect under which the aerodynamic coefficients are not 
uniquely dependent upon the instantaneous incidence angle but rather depend 
upon the time history of the airfoil motion. These differences between static 
and dynamic stall indicate that even though Btatic and dynamic stall may have 
similar mechanisms, theories which employ static stall data to predict dynamic 
stall phenomena cannot be expected to be accurate. Accurate theoretical 
predictions of the airfoil loading during dynamic stall require a theory which 
recognizes the time -dependent nature of the dynamic stall problem. 
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Several theories of various degrees of sophistication have been proposed 
to predict airfoil loading during dynamic stall. One such theory, due to Ham 
(ref. 4), is based upon an inviscid flow analysis which ignores any direct 
effect of boundary layer separation upon stall phenomena. The theory models 
dynamic stall by the shedding of vortices from the airfoil leading edge. 
Although the theory has predicted both the maximum lift and moment coeffi- 
cients for an airfoil undergoing dynamic stall, it has not yet predicted the 
coefficients through an entire loop. In addition to Ham's inviscid theory, 
semiempirical dynamic stall theories are also available. A semiempirical 
method used by Carta, Commerford, and Carlson (ref. 5) is based upon a 
correlation of existing experimental data. However, it is not clear how far 
the method can be extended to either other airfoils or other types of motion. 
Similar procedures could be developed for airfoils undergoing other types of 
motion, however, each class of airfoil pnd motion may require a different body 
of experimental data. The procedure of Ericsson and Reding (refs. 6 and 7) is 
based upon assuming an effective camber and effective incidence which 
change as the hysteresis loop develops. When these effective quantities are 
used in conjunction with a semiempirical lag time, the force and moment 
coefficients during dynamic stall are predicted. However, the method is 
highly-empirical and does not predict stall from basic boundary layer 
separation considerations . 

In contrast to the analyses of refs. 4 through 7, which are semiempirical, 
there is the more fundamental analysis developed by Crimi and Reeves (ref. 8). 
Hie Crimi -Reeves analysis is based upon a solution of the fluid dynamic 
equations in the neighborhood of an airfoil in arbitrary motion. In brief, 
the analysis of ref. 8 combines the solution of the linearized potential flow 
equations with the boundary layer momentum equations to predict the flow field 
behavior . Although the procedure produces qualitative agreement with the 
basic features of dynamic stall, its theoretical predictions are in quantita- 
tive disagreement with experimental data. An exrmination of the Crimi -Reeves 
analysis indicates assumptions are made that may lead to the observed 
quantitative differences between theory and experiment. In particular, the 
Crimi -Reeves analysis is based upon simplified treatments of separated regions, 
the transition process, and the nominally inviscid flow field. Although the 
method of ref. 8 uses a finite -difference solution to the boundary layer 
equations in regions of attached flow, it uses an integral boundary layer 
solution in regions of separated flow. The integral solution requires an 
assumption of a velocity profile family and this assumption restricts the 
validity of the solution. Secondly, the procedure uses an empirical transition 
model and, finally, the procedure uses a linearized potential flow theory to 
represent the outer inviscid flow. 
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The present report describes the development of a weak-interaction 
solution for the dynamic stall of helicopter rotor blades using an approach 
similar to that of Crimi and Reeves (ref. 8), but with an improved treatment 
of the separated flow regions, transition phenomena, and the potential flow 
regions. By definition, the weak-interaction solution ignores the effect of 
viscous displacement upon the nominally inviscid outer solution. The weak- 
interaction assumption is valid as long as the viscous displacement thickness 
remains small compared to the airfoil thickness. However, when the displace- 
ment thickness becomes large and significantly modifies the inviscid pressure 
distribution, such as in regions of significant boundary layer separation, the 
theory is invalid. In order to obtain accurate predictions of the flow field 
when a significant separation region is present, it is necessary to use a 
strong -interaction theory which recognizes the mutual interaction between the 
viscous inner and nominally inviscid outer flow fields. Such a strong- 
interaction calculation procedure could be developed by an extension of a suc- 
cessful weak-interaction procedure in which an inner solution such as the 
the viscous solution of the present report is coupled to an inviscid outer 
solution. The coupling would require continuity of flow angle along the line 
joining these solutions. Alternatively, the entire flow field could be 
solved by the Navier -Stokes equations thus avoiding the problem of coupling 
two different solutions in two regions of the flow. Such a solution has 
recently been obtained for internal duct flow problems by Briley and McDonald 
(ref. 9). 

Although the present weak-interaction solution is limited in 
applicability to flow situations in which the viscous displacement thickness 
does not significantly afxect the inviscid pressure distribution, the present 
effort can accurately predict the flow field under conditions for which the 
boundary layer does not significantly affect the pressure distribution. In 
this regard, as is shown subsequently, the procedure is capable of analyzing 
the detailed viscous flow mechanisms including those mechanisms governing the 
leading edge separation bubble and the method is also capable of predicting 
conditions for incipient stall. In addition, the procedure may be regarded as 
a first step in the development of strong-interaction solution to the isolated 
airfoil dynamic stall problem, or a solution for the dynamic stall problem 
based upon a single set of equations representing the flow field in the entire 
solution domain. 

Hie authors are pleased to acknowledge the considerable assistance 
contributed by Dr. W. R. Briley of United Aircraft Research Laboratories to 
this effort. Dr. Briley furnished the authors with a detailed explanation of 
the viscous flow computer code and contributed to many of the ideas presented 
in this report through on-going discussion of numerical calculations of 
viscous flow fields. 
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LIST OF SYMBOLS 

structural coefficients 
chord 

skin friction coefficient 

sublayer damping function 

low Reynolds number correction function 

turbulence source terms 

curvature, or reduced frequency 

mixing length 

dissipation length 

pressure 

turbulence kinetic energy 
turbulence Reynolds number 
turbulence source terms 
surface coordinate 

surface coordinate of stagnation point 
time 

streamwise velocity component 
transverse velocity component 
streamwise coordinate, or chordwise location 
transverse coordinate 


incidence angle 


boundary layer thickness, or calculation layer thickness 


sublayer thickness 


displacement thickness 
turbulence dissipation 


dimensionless transverse coordinate 


momentum thickness 


von-Karman constant 


kinematic viscosity 


kinematic eddy viscosity 
vorticity 


density 


acceleration parameter 


shear stress or relaxation time 


integral functions (see Eqs. (35) through (38)) 


stream function 
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THEORY 


General 

Hie airfoil flow field calculations described in the present , i h • 
obtained by dividing the flow field into several distinct regions c • solvii. . 
the appropriate equations in each region. In this manner results obtained in 
one region serve as boundary conditions for other regions and the entire flow 
field about the airfoil is constructed. The major flow field division i l- the 
separation between the relatively th -i viscous layer near the airfoil surface 
and the nominally inviscid outer flow. The calculations are initiated by 
obtaining the nonlinear, incompressible, inviscid solution about an airfoil 
in arbitrary motion using the procedure of Giesing (ref. 10). The Giesing 
procedure predicts an inviscid velocity distribution about the airfoil as a 
function of time and this inviscid velocity distribution serves as a time- 
dependent outer edge boundary condition for the viscous flow calculation. The 
viscous flow region is divided into several subregions, as shown in Fig. 1. 

Hie stagnation region is specified in the vicinity of the airfoil leading edge. 
This is followed on the si tion surface by the region where a leading edge 
separation bubble is anticipated and by a fully-turbulent region which may 
exhibit a trailing edge separation bubble. The fully-turbulent region may, 
for convenience, be divided into two or more subregions. The stagnation 
region is followed on the pressure side by a transition region where boundary 
layer transition is likely to occur and then by a fully-turbulent region which 
again may be divided into two or more subregions. 

The viscous calculation is initiated in the stagnation region using the 
inviscid flow solution as an outer edge boundary condition. The necessary 
boundary conditions at the junction between the stagnation and transition 
region and the stagnation and possible separated region are specified by 
assuming that at these boundaries the second derivatives in the streamwise 
direction are zero. When the viscous flow field is divided into several 
distinct subregions, a significant advantage is gained in terms of the 
required computer storage. Obviously, for a given grid resolution the storage 
required by the subregion approach is much less than that which would be 
required if the entire viscous flow were done in a single calculation. However, 
if a large enough core capacity were available, the subregion approach would 
not be necessary. The subregion approach has the disadvantage of limiting the 
amoun , of upstream Influence, Since the upstream boundary conditions for each 
region (except for the stagnation region) are determined by the flow field of 
the previous region, it is obvious that upstream influence cannot propagate 
through a subregion boundary. This constraint upon upstream influence can 
lead to a constraint on the upstream propagation of a separated flow region as 
demonstrated in the airfoil calculation discussed subsequently. 
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After a solution is obtained for the stagnation region, solutions are 
found for the separated and transition regions. In each case the upstream 
boundary condition is obtained from the stagnation region solution; at the 
upstream boundary vorticity and stream function are continuous. The outer 
edge boundary condition is imposed by the inviscid flow solution and the 
downstream boundary condition is set using the assumption that second 
derivatives with respect to the streamwise coordinate are zero. The solutions 
in the fully-turbulent region are then obtained in a similar manner. 


The Potential Flow Solution 

The outer nominally inviscid, potential flow field is obtained using the 
computer code developed by Giesing (ref. 10). The Giesing procedure calculates 
the incompressible, inviscid flow field about a two-dimensional airfoil in 
arbitrary unsteady motion. The basic technique used by the procedure is to 
divide the airfoil body surface into a series of elements and then to apply a 
source distribution around the body, one at each element. The source distri- 
bution is adjusted until the velocity normal to the airfoil surface is zero. 

In addition to the airfoil source distribution, the airfoil contains a bound 
circulation which is adjusted to satisfy the Kutta condition at the airfoil 
trailing edge. The vortex wake shed by the airfoil is carried by the fluid 
particles to which it is attached and changes i*' space and time. A complete, 
detailed description of the procedure and coopai sons with experimental data 
and other analyses can be found in refs. 10 and xl. 


The Viscous Flow Solution 

General . - The viscous calculation procedure used in the present report 
was originally developed by Briley and McDonald (ref. 12) in a study of 
transitional separation bubbles. The procedure solves the vorticity-stream 
function incompressible Navier-Stokes equations in either their full or 
reduced form. When the full vorticity-stream function form of the equations 
is used, the time -dependent vorticity equation is solved with a Douglas-Gunn 
(ref. 13) perturbation of the Crank -Nicholson differencing scheme which 
generates an alternating-direction -implicit (ADI) method. In the reduced form, 
a Douglas-Gunn (ref. 13) perturbation of the backward difference scheme is 
used for the vorticity transport equation while an approximate stream function 
equation, not requiring a full ADI sweep, is used. In both the full and 
approximate form of the governing equations the resulting ADI procedure 
advances the vorticity equation in time through a two-step calculation 
procedure. ADI procedures present an extremely efficient way of solving multi- 
dimensional problems, and in the past a major obstacle in the routine solution 
of such problems has been the large amount of cos^uter time required to obtain 
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a solution. In explicit methods such as those of ref. 14 the governing 
equations are subject to one or more stability restrictions on tS ■? size of 
the time step relative to the spatial mesh size. These stability limits 
usually correspond to a limit arising from convective considerations or to a 
limit arising from diffusive considerations. In typical incompressible 
boundary layer type flows the viscous limit is expected to dominate the 
problem. 

While explicit stability restrictions may not be a serious disadvantage 
for inviscid flow calculations in situations for which a laminar boundary 
layer must be resolved the limitation in time step specified by the explicit 
procedure stability limits may be a serious disadvantage. In the case of a 
turbulent boundary layer the extreme resolution required in the sublayer 
region obviously magnifies the problem. Thus, explicit methods inherently 
contain the key disadvantage that the maximum step size is fixed by the 
spatial step size rather than the physical time scale. One way ort of the 
sublayer dilema with an explicit procedure is simply to assume a form for the 
appropriate variable within the sublayer. In the present problem of flow 
separation it is obviously very limiting to suggest such a profile as an 
artifice to eliminate grid points, particularly in a flow where the temporal 
behavior close to separation is of interest. 

In contrast to explicit methods, implicit methods tend to be stable for 
large time steps anc so have the allowable time step set by physical consider- 
ations rather than by the computational mesh and, therefore, offer the 
prospect of substantial increases in computational efficiency if the computa- 
tional effort per time step is competitive with that of explicit methods. In 
the calculations of the present report the viscous stability limit would be 
expected to determine the maximum time step limitation of an explicit 
solution method, hcwever, due to the use of an implicit method, time steps on 
the order of 1000 times the viscous stability limit were not uncommon. 

The governing equations. - The calculations presented in tr.e present 
report are based upon a solution of the Navier -Stokes equations written in 
vorticity-stream function form. The equations are written in an airfoil 
coordinate system for airfoil sections subjected to two-dimensional flows. 
Within these two-dimensional boundary layers Coriolis and centripetal effects 
due to the pitching motion of the airfoil are expected to be small compared to 
the viscous and usual boundary layer convective effects and, therefore, they 
are neglected. It should be noted that at the leading edge stagnation point 
Coriolis effects may be important, h"**ver, these effects became negligible 
very quickly. For example, if the leading edge region of the NAC k 0011 air- 
foil is approximated by a circular cylinder, a : faced frequency of 0,2 and a 
chord Reynolds number of 10? is assumed, and for the purpose of estimation of 
transverse velocity gradients, a steady flew calculation is used to approlmate 
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the viscous flow near the front stagnation point then the Coriolis force 
becomes negligible compared to the convective momentum transfer well within 
one degree of the stagnation point along the cylinder surface. As shown by 
many authors (e.g., ref. 15), the vorticity transport equation can be written 
in Cartesian coordinates as 


at 3 % 


+ 


dy 




( 1 ) 


where | is the vorticity given by 


du dv 
dy dx 


( 2 ) 


t is time, x and y are the Cartesian coordinates, u and v are velocity 
components in the x and y directions, respectively, and v is the kinematic 
viscosity. A stream function is defined by 


, — 

dy 


v 


d* 

dx 


( 3 ) 


which leads to the relation 
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Equations (l) and ( 4 ) form the set of equations to be solved for laminar 
flow. For turbulent flew the equations are written in the form of conservation 
of momentum in the x- and y-coordinate directions as 


du du du 
— + u — tv — 

dt dx dy 


i it + Ji!i + iii] 

P dx Idx 2 dy 2 J 


( 5 ) 
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( 6 ) 


dv dv dv 

— + U — + V 

dt dx dy 


iE + 4A 

dy Idx 2 



The velocity components u and v are divided into mean and fluctuating parts 


u * 0 + u 


(7) 


v 


V + v' 


( 8 ) ' 


where the overbar indicates mean quantities and the prime indicates fluctuating 
quantities. Substituting Eqs. (7) and (8) into Eqs. (5) and (6) and averaging 
leads to 
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where 
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The stress -strain relations are assumed to be given by 


r<r x r xy | |o(u y + v x ) u y + v x | 

lr xy <r y j P f l u y ♦ v x b(u y ♦ v x )J 


(12) 


where is the kinematic eddy viscosity and the subscripts x and y indicate 
partial differentiation with respect to x and y, respectively. When Eq. (12) 
is substituted into Eqs. (9) and (10), use is made of the relation 


c m 6± .*1 ( 13 ) 

dy d% 

and the pressure is eliminated from the equations, the resulting turbulent 
vorticity transport equation is 



where S is the turbulent dissipative contribution given by 

s - • 3f„ - 47,,,] + (f ♦ 25.) 
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(HO 


(15) 


In the calculations made in the present report, the viscous flow fields are 
basically boundary layer type flews in which the length scale and velocity in 
the streamwise direction are at least an order of magnitude larger than the 
length scale and velocity in the transverse direction. If the x-coordinate is 
associated with the streamwise direction, this implies that 
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so that 


= *T fyy + *tyy ^ * 2 *ty ^y 


(17) 


For turbulent flow Eqs. (4), (l4), and (17) determine the flow field. 

Equations (l) and (4) for laminar flow or Eqs. (4), (l4), and (17) for 
turbulent flow represent the full Navier -Stokes equations (with the approxima- 
tion to turbulence dissipation) in stream function-vorticity form. If the 
flow field being investigated is a boundary layer type flow in which the 
streamwise velocity, u, is much larger than the transverse velocity, v, and 
derivatives with respect to the transverse coordinate, y, are much larger than 
derivatives with respect to the streetwise coordinate, x, then 



du 

*y 


(18) 


and the stream function equation, Eq. (4), may be approximated by 



(19) 


For these flows Eq. (19) replaces Eq. (4). When Eq. (19) replaces Eq. (4), 
the set of equations is termed the reduced set of equations. It should be 
noted that the reduced set of equations is equivalent to a set of boundary 
layer equations with the addition of a streamwise diffusion of vorticity term 
(ref. 12). In most of the calculations presented in the present report the 
streamwise diffusion term does not contribute significantly to the vorticity 
equation balance. 


l4 



Hie preceding equations represent the equations written in a Cartesian 
system; however, for purposes of the calculations in the present, report, a 
coordinate transformation was performed to allow for airfoil curvature and so 
as to better represent the change in shear layer width along the streamwise 
direction. The curvature of the body, k, is accounted for by assuming a 
coordinate system for curved walls in which the streamwise distance, x, is 
taken along the body surface and the transverse distance, y, is taken normal 
to the body surface. The curvature of the body, k, is a specified function of 
x. The curved-wall coordinate system is discussed in detail in ref. 15 . In 
the curved -wall body coordinate system the vcrticity equation becomes 



where 5 is vorticity, u and v are velocities in the x and y directions, 
respectively, and v is kinematic viscosity. The vorticity-stream function 
relation remains 


V 2 * * t 


( 21 ) 


which in the curved-wall coordinate system is written as 


I dfy 

( l + ky ) 2 <** 2 


y dk + d 2 \jf + k 

(1+ky) 3 dx dx dy 2 (i + ky) 


s *( 22) 

dy 


Tht velocities are related to the stream function by 


U* 



( 23 ) 
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V* 


(24) 


- 1 dty 
i+ky dx 


For turbulent flew the momentum conservation equations in the x- and y- 
directions are written in the curvilinear coordinate system and the velocity 
is divided into mean and fluctuating parts. After the usual averaging 
procedure is performed, the pressure is eliminated by cross differentiation 
and an order of magnitude argument is applied to the apparent shear stress. 
The procedure is analogous to that of 
result 

af u af . af 

at + l ♦ ky dx + V ay 

v dk d£ 

( i + ky) s y dx dx 


The derivation of Eq. (25) is given in APPENDIX A. Equations (22) and (25) 
represent the Navier-Stokes equations in curvilinear coordinates. In order to 
always keep the grid within the viscous layer a further transformation is 
introduced whereby 


Eqs. (5) through (15), and leads to the 


v d z i a 2 
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(25) 


(l + ky) ay 


x = x 


(26) 


V 3 y/8 


(27) 
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The vorticity transport equation becomes 





where subscripts indicate differentiation and 


Z *1+ krjB (29) 

The stream function equation becomes 


It should be noted that the final form of the equations allows the grid to 
adjust naturally to both spatial and temporal changes in the boundary layer 
thickness. 
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Hie Turbulence and Transition Models 


The turbulence model. - In laminar flows the stress is composed solely of 
the molecular stress which is determined uniquely by the molecular viscosity 
(a property of the fluid) and the local velocity field. However, as shown in 
Eqs. (9) and (10), in turbulent or transitional flow where the flow field is 
composed of a mean and fluctuating part the averaging procedure gives rise to 
additional terms which appear to act as stress terms in the momentum conserva- 
tion equations and which generally are called turbulent stress terms. These 
additional terms in the momentum equation lead to additional terms in the 
vorticity transport equation and calculations in the turbulent and 
transitional regime require a mathematical model to represent these terms. 

Since, under normal conditions, airfoils operate in the regime where laminar, 
transitional, and turbulent viscous flow regions are present, any calculation 
procedure attempting to predict the viscous flow field about an airfoil must 
contain both a turbulence model and a transition model. 

Insofar as the turbulent flow is concerned a large variety of models have 
been developed for fully-turbulent flows (e.g., ref. ^5). These models can be 
divided into two broad categories, equilibrium turbulence models and historical 
turbulence models. Equilibrium turbulence models assume that the turbulent 
stress is determined uniquely ty the local mean velocity field. These equilib- 
rium procedures usually hypothesise an eddy viscosity or mixing length deter- 
mined solely be mean flow conditions. Although equilibrium turbulence models are 
adequate for the prediction of many turbulent boundary layers, their basic 
assumption relating the turbulent shear stress to local flow conditions is 
clc \rly in error for flow situations in which rapid changes in the mean flew 
field occur (e.g., refs. 17 and l8); in rapidly developing flows the turbulent 
stress is not determined by local conditions but rather by the history (both 
upstream spatial history and temporal history) of the flow and in these cases 
the theory is improved if a model which includes the flew history is used. A 
large number of such models for steady-state turbulent flows are discussed in 
ref. 16. In addition to the steady-state historical turbulence models of 
ref. 16, a time -dependent historical model has been developed by Patel and 
Nash (ref. 19) and applied to a variety of flows by Nash, Can, and Singleton 
(ref. 20). 

In the present calculations the turbulence model used is that developed 
by McDonald and Camarata (ref. 21) which solves an integral form of the 
turbulence kinetic energy equation. In brief, the turbulence kinetic energy 
equation is a conservation equation derived from the Navier-Stokes equations 
by writing the instantaneous quantities as a sum of mean and fluctuating parts. 
The ith Navier-Stokes momentum conservation equation (i = 1,2,3 referring to 
the three coordinate directions) is multiplied by the i^ component of 
fluctuating velocity and the average of the resulting three equations is taken. 
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The three averaged equations are summed to obtain the turbulence kinetic energy 
equation. The derivation of the turbulence kinetic energy equation has been 
given by Favre (ref. 22) for the general case of a compressible fluid and 
approximated by Bradshaw and Ferris (ref. 22 ) to boundary layer flows. For 
incompressible flow the derivation has been given by many authors (e.g., ref. 
24). 


In the case of incompressible flow, the boundary layer approximation to 
the turbulent kinetic energy equation can be written as 


■jr (**A) ♦ +-£(v<M) « -77 -4L 

temporal changs advsction production 

(31) 

diffusion dissipation normal itrsss production 

where q 2 is the turbulence kinetic energy. 



■ 7 * ♦ 7 * ♦ 



(32) 


It should be noted that in the turbulent flow calculations of the present 
report the boundary layer approximations are expected to be valid and thus 
the boundary layer form of the turbulence kinetic energy equation is 
appropriate . 

Following Townsend (ref. 25) and Bradshaw and Ferris (ref. 23), McDonald 
and Camarata (ref. 21) introduce structural coefficients a n and L, together 
with a mixing length £; chese scales are defined as 


u' v' 


« 


■ o,y , 7 1 ■ a,?" , 7* *0,? , 

.C77J * / 7l 


( 33 ) 


19 



and a^, & 2 > an ^ a 3 ar * assumed to have values 0.15, 0.50, and 0.20, 
respectively. The turbulence kinetic energy equation is then integrated 
between the wall and the edge of the boundary layer at y = 6 to obtain the 
equation 

♦ E (34 > 

where 

*'‘+C* (35) 

**’ t C x s *(■$ *> (36) 

+a * * [• - * L ]|*J'|* *» (37) 

+* * 2 *' |-Jj- 1 [iF *» #] 6r > (38) 




( 39 ) 



\ 


( 40 ) 

l 

Hie term containing 0i represents the temporal rate of change of turbulence 
- kinetic energy, the term containing represents the streamwise rate of 

change of turbulence kinetic energy, the term containing 0^ represents the 
integral of turbulence production minus dissipation, and the term containing 
04 represents the normal stress production. The terms designated by E 
represent turbulent source terms resulting from disturbances imposed upon the 
viscous layer by the free stream. It should be noted that the turbulence 
kinetic energy at the edge of the boundary layer, q2 e is damped by the factor 
U e A®) 2 where l e is the value of the mixing length at the boundary layer edge 
and I,,, is the ,T wake" value of mixing length; i.e., the value far from the 
wall. In most regions of the flow ljl m =1.0 and thus no damping of the 
. entrained turbulence energy occurs. However, near the stagnation point, £ e 

is initially considerably smaller than j l* due to the extremely thin highly 
viscous layer present near the stagnation point and the entrainment is damped 
„ heavily in this region. As the flow proceeds away from the stagnation point 

r le quickly rises to the value of and no damping of the entrained turbulence 

occurs. 
r 

^ Hie dissipation length is given by 

[ l = oiB tonh [«y / (o.l 8 >] (4i) 


I 


I 

* 

1 


where K is the von Karman constant taken as 0.4l, is a sublayer damping 

factor, and is a low Reynolds number correction. In the original McDonald - 

Camarata model the sublayer damping was assumed to be given by the van Driest 
damping model and no low Reynolds number correction was made. However, 
following a later work by McDonald and Fish (ref. 26), the su blaye r damping is 
assumed to distribute normally about a mean height y + (y + = y^T/p/v) with a 
standard deviation o leading to the equation 

(y + -y*) /«■} <■*> 


where P is the normal probability function; y* is taken as 23, and a as 8. 



The low Reynolds number correction is based upon the work of McDonald 
(ref. 27) which relates the correction factor, 2b 2 » to a turbulence Reynolds 
number, R T , given by 


R r » 



(43) 


where 6 is the viscous layer thickness, 6 S is the sublayer thickness defined 
as the location at which the laminar stress has fallen to 4 percent of the 
total stress, v^. is the turbulent kinematic viscosity defined as 

■ ( “ul v')/(dO/dy) (44) 


and v is the laminar kinematic viscosity. The correction factor ^ 2 is given 
by 


2> t « [l.O + txp ( - 1.63 In f T + 9.7)] 


(45) 


where 


f T * 681 R t + 614.3 R r >40 


(46) 


f T «ioo* T ft r t si 


(47) 


and for l<H T s4o, f T is fitted between Eqs. (46) and (47) by a cubic 
constructed to match the function and slope at the join points. 
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Finally, a one-parameter mixing length profile, 4* , is introduced where 


/ • ^ 2, txn [ify/ij 


( 48 ) 


Introduction of Eqs. (4l) and (48) into Eq. (34) yields a differential equation 
for the wake value of the mixing length, i m , which is solved in conjunction 
with the mean flow equations to determine both the mean flow field and the 
shear stress development. A detailed derivation and discussion of the turbu- 
lence model as well as the transition model (to be discussed subsequently) is 
given by Shamroth and McDonald (ref. 28). 

The transition model. - Althougn the McDonald -Camara- " -bulence model, 
as previously described, is a well -proven turbulence model for fully-turbulent 
flows, it is still necessary to include a model to predict the flows in the 
transitional regime. Such a model based upon a solution of the turbulence 
kinetic energy equation has been developed by McDonald and Fish (ref. 26) and 
has been verified through comparisons with a large body of experimental data 
by McDonald and Fish (ref. 26), Shamroth and McDonald (ref. 28), and 
Kreskovsky, Shamroth, and McDonald (ref. 29). It should be noted that this 
transition model is based upon a rigorous conservation equation rather than 
semiempirical data correlations, as is the case with most other transition 
models (see, ref. 30). The model has successfully predicted the behavior of 
a large variety of transitional boundary layers from the incompressible to the 
low hypersonic Mach number regime subject to various heat transfer rates, 
pressure gradients, and wall roughness heights (refs. 26, 28, and 29). In 
brief, the transition model is Identical to the turbulence model with the 
exception that the structural coefficient, a^, becomes a function of the 
turbulence Reynolds number of the form 

V°o [ , t'\]/|'+«**<\> [(» r /f To ) -.]j (49) 

where f T is given by Eqs. (46), (47), and the cubic fit with f To = 100. The 
variable ia a function of the wall-to-free-stream temperature ratio (ref. 

28) and for the case of wall temperature equal to free stream temperature a 0 
is equal to 0.011$. Thus the turbulence kinetic energy equation is always 
solved in conjunction with the governing mean flow equations. If the calcu- 
lated mixing length is very small, the flew is laminar; however, the mixing 
length may increase causing the turbulent transport to be comparable to the 
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lamina:r transport and, in this case, Lhe flow is transitional. Finally, the 
mixing length may reach a fully -turbulent value leading to a fully-turbulent 
flow. 


Solution of the Equations 

Solution of the full set of equations. - When the full set of vorticity- 
stream function equations are used, an iterative solution is required between 
the vorticity transport equation, Eq. (l), and the vorticity-streej function 
relation, Eq. (4). Hie set of equations are subject to the boundary condi- 
tions at the wall and at the boundary layer outer edge as follows ; 


at y = 0 


<?- 

n 

0 

(50) 

30/a y = 0 


at y = 6 


30/3 y = Ug 

(51) 


§ = 0 


Briley and McDonald (ref. 12) solve the set by first assuming a wall vorticity 
distribution and then solving the vorticity transport equation by a Douglas - 
Gunn ADI procedure. In the Douglas -Gunn procedure the vorticity transport 
equation is written in the general form 


*•■**„♦ +C V°* 



(52) 


where 4 is a dependent variable (in this case, vorticity), A, B, C, D, and E 
are coefficients, and subscripts indicate derivatives. A two-step calculation 
procedure is med to advance the calculation in time from t 0 to t 0 + At. In 
this procedure 3 the dependent variable at the town time step, t 0 , 4 >* i^ 
the dependent variable at the first step of the calculation procedure, and 
0** is the dependent variable at the second step of the calculation (or at 
time = t 0 + At). 
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In the first step of the procedure derivatives with respect to y are 
treated implicitly using Crank-Nicholson differencing and those with respect 
to x are treated explicitly leading to the equation 

Af ' ■§■ & 4 4 4 £ 4 *f ♦» 4 E (53) 


which is solved for <t>* by inverting a tridiagonal matrix. In the second step 
of the procedure derivatives with respect to x are treated implicitly while 
those with respect to y are treated explicitly leading to the equation 


£--£■ *7 - £ £ +‘„- i*‘»= -fr 

4 "S’ ♦m 4 + 4 T*y 4E 


Equation (53) is subtracted from Eq. (5*0 leading to the simpler equation for 
</>** 




(55) 


After the vorticity transport equation is solved, the stream function 
equation is solved with the wall and edge boundary conditions 


0 = 0 at y = 0 
30/3y = u e at y = 6 


(56) 


by a Peaceman-Rachford ADI procedure for elliptic equations (ref. 31). The 
Peaceman-Rachford procedure introduces a series of acceleration parameters, p^, 
and then solves the stream function equation by an iteration process. For each 
choice of acceleration parameter, pi, two sweeps are made. In the first sweep 
the iteration is advanced from solution to by treating x-derivatives 
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as explicit and y-derivatives as implicit. In the second sweep the x -derive - 
tives are treated implicitly and the y-derivatives explicitly and the solution 
is advanced to Thus, the procedure advances as 


i+l/2 

yy 


♦ ft,* 


i ♦ •/* 




,i+i ,i+i 

*«« ♦ *+ • 


* y yy 




1 + 1/2 


(57) 


If agrees with i// 1 within a specified tolerance, the equation is 

considered solved; if not, the procedure is repeated for the next acceleration 
parameter. It should be noted that, if the acceleration parameter, p^, is 
taken as the inverse of a time step l/At^, the Peaceman-Rachford procedure is 
analogous to the Douglas -Gunn procedure. 


After the stream function equation has been converged the streamwise 
velocity component at the wall consistent with the stream function solution is 
computed. If the wall velocity is zero within a specified tolerance, the set 
of equations is considered solved at the given time step; if not, a new wall 
vorticity distribution is assumed using a secant extrapolation and the entire 
procedure consisting of the solution of the vorticity transport equation and 
the stream function -vorticity relation is repeated. 

Solution of the reduced set of equations. - As can be summarized from 
the preceding discussion, the iterative solution to the full set of equations 
can be a relatively time consuming process. The computer time required by the 
reduced equations is less by an order of magnitude than that required by the 
full set of equations; therefore, the reduced set represented by Eqs. (l) and 
(19) have been utilized whenever possible. When the reduced set of equations 
are used, the vorticity transport equation and the reduced stream function 
equation are solved as a coupled set by the Douglas -Gunn ADI perturbation of 
the backward difference procedure using the boundary conditions represented 
by Eqs. (50) and (51). When the reduced set of equations are solved, the 
equations are first integrated with the x -derivatives taken as implicit and the 
y-derivatives as explicit. Since the reduced stream function-vorticity rela- 
tion contains no derivatives with respect to x, the integration is performed 
only for the vorticity transport equation and is similar to that for the full 
set of equations. The equations are then integrated with y-derivatives taken 
as implicit and x-derivatives as explicit. During these integrations a 
coupled set of equations is solved in which the coefficient matrix takes the 
form of a block tridiagonal matrix, each block being a 2 x 2 submatrix. The 
fact that the system can now be solved in a coupled manner eliminates the 
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previously required iteration on wall vorticity. The reduction in the stream 
function equation eliminates the Peaceman-Rachford iterative ADI solution 
required when using the full stream function equation. 


RESULTS 


Verification of the Calculation Procedure 

General . - Before predictions were made for the flow field about 
oscillating isolated airfoils, it was deeded necessary to verify the calcula- 
tion procedure. Since the potential flow computer code used in the present 
effort was extensively checked out by its originator (refs. 10 and 11), little 
effort was expended in verifying this code. The major verification effort 
concentrated on the viscous flow calculation procedure. Verification of the 
viscous flow calculation procedure ,'equires verification of the mathematical 
models used in the procedure and verification of the finite -difference proce- 
dure itself. Hie finite -difference procedure was verified by comparing 
predicted results and available analytical solutions for a wide variety of 
laminar flows. These comparisons between numerical solutions of the present 
procedure and analytical solutions for laminar flows assess how well the 
finite -difference solutions correspond to solutions of the original differ- 
ential equations. In addition, calculations for transitional boundary layers 
made with the present procedure were compared with calculations made with the 
well-established UARL boundary layer prediction deck. 

Insofar as the mathematical models are concerned, the only mathematical 
models present other than that implied by representing the flow field by the 
Navier-Stokes equations, are the turbulence and transition models. These 
models have been verified for steady-state, unseparated boundary layer type 
flows through a large number of comparisons between theoretical predictions 
and experimental data (refs. 26, 28, 29, and 32). In addition, successful 
predictions of experimentally measured transitional separation bubbles have 
been made by Briley and McDonald (ref. 12). Thus, the turbulence and transi- 
tion models are well-established for steady-state viscous flows. 

The major unresolved question centers upon how well the turbulence and 
transition models represent time -dependent flew fields. Since the basic 
equation used to predict the turbulent shear stress development is the 
turbulent kinetic energy equation, an equation derived directly from the 
Navier-Stokes equations, there can be little argument regarding the governing 
equation itself, Eq. (31) . The major uncertainty resides in the turbulence 
model, Eqs. (33) and (4l) through (48); however, some such model nrst be 
hypothesized in order to solve the turbulence kinetic energy equation. If 
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unsteady calculations are to be valid, then the turbulence model should be 
valid for the flows under consideration. In general, the characteristic 
frequency of turbulence is given by 




CjQ) 


and for helicopter applications the characteristic frequency of the airfoil is 
at most 


<*»<,* 


* 
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where 6 is the boundary layer thickness, u is the free stream velocity, and c 
is the airfoil chord. The ratio of the frequency of the turbulence to that of 
the airfoil motion is 



= JL 

28 


>>l 


(60) 


Therefore, the turbulence frequency is expected to be much greater than the 
airfoil frequency and the turbulence structure should be unaffected by the 
time -dependent mean motion at least in the absence of large separated regions. 
Thus, the turbulence model is expected to be valid. 

Potential flew calculation. - As previously discussed, the potential flew 
calculation procedure is a well accepted procedure, the results of which have 
been verified by its author (refs. 10 and 11). An additional comparison was 
made under the present effort between theoretically predicted lift and moment 
coefficients with the data of Carta, Commerford, Carlson, and Blackwell (ref. 
33). The results of this comparison, shown in Fig. 2, are considered to be 
good. Hie experimental data have been corrected for finite span effects by a 
correction factor obtained from comparison of the low Incidence experimentally 
determined lift -incidence slope with the theoretical value. 
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Laminar viscous flow calculations. - One suitable test case for assessing 
the present calculation procedure is the flat plate stagnation point laminar 
flow field termed Hiemenz flow (ref. 15). The inviscid flow solution at large 
distances from the plate is given by 


u = ax 
v = -ay 


(61) 


The analytical viscous solution is a similarity solution in wh’ rh t v '-' 
dimensionless velocity, u/ >, is a function of the dimensionless transverse 
coordinate, y/./v/a , where Ug is the inviscid flow velocity in the streamwise 
direction, and v is the kinematic viscosity. A comparison ottween the 
solution predicted by the numerical proceduit and the analytical solution as 
given by ref. 15 is presented in Fig. 3. As can be seen, the comparison is 
excellent . 

A second steady-state laminar solution which serves as a good test case 
for the calculation procedure is that of flow about a circular cylinder, since 
such a flow provides an additional feature not present in the Hiemenz flow. 

In particular, flow about a circular cylinder provides a geometry having a 
finite -radius of curvature and thus the calculation serves as a check on the 
radius of curvature effect in the present calculation procedure. A comparison 
between the numerical solution of the present analysis md that obtained 
analytically from the Blasius series solution to the boundary layer equations 
(ref. 15) is presented in Fig. 4. The agreement between the numerical and 
analytical solutions, shown in !*lg. 4, at three different angular locations 
(8=0 being the front stagnation point) is excellent. In Fig. 4, ue is the 
local free-stream velocity, is the approach velocity, and R is the 
cylinder radius. 

The viscous flow calculation procedure was also used to predict the flow 
along a flat plate oscillating sinusoidally in its own plane to assess the 
accuracy of the procedure in predicting time -dependent flows. The numerical 
solution was run using both the full set of Navier -Stokes equations and the 
reduced set of equations; no significant differences in the results were 
apparent. Hie velocity profiles calculated by the numerical procedure through 
one half a cycle are compared to the analytical solution of ref. 15 in Fig. 5. 
Hie analytical solution is given by 
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u(y . t ) - u 0 [cos M) - e* hy cos (ut - ky) 


( 62 ) 



(63) 


In Fig. 5 the open circles represent the solutions obtained from both the 
reduced and full Navier-Stokes equations for an oscillation frequency of 
100 radians/sec. It should be noted that in the numerical calculation the 
plate was assumed stationary and the inviscid flow field oscillated, whereas 
in ref. 15 the plate was assumed to oscillate in a quiescent fluid; however, 
as discussed by Lighthill (ref. 34), the two problems are equivalent. As can 
be seen, agreement is excellent. A comparison between the wall vorticity, 
3u/3y, predicted by the numerical procedure and that given by the analytical 
solution is presented in Fig. 6. The wall vorticity, being proportional to 
the skin friction, is a very sensitive indicator of the accuracy of the 
numerical procedure. The agreement is again excellent. It should be noted 
that the reduced equations were used successfully for the oscillating plate 
problem even though regions of reversed flow are present in the flow field. 

In addition to these laminar viscous flow calculations made in the 
presence of a solid wall, calculations were made for the wake behind a flat 
plate and these results are compared to the results obtained by Goldstein and 
Luckert as given by Rosenhead and Simpson (ref. 35) in Fig. 7. Since 
Goldstein's solution was numerical and Luckert 's graphical, Goldstein's 
solution is regarded as the more reliable of the two. Agreement between the 
present solution and that of Goldstein is good. 

Transitional viscous flow calculations. - The results of Figs. 3 through 
7 demonstrate the ability of the present numerical procedure to predict 
boundary layer, stagnation point, and wake flows in the laminar regime 
including unsteady effects and the effects of a finite radius of curvature. 
The procedure's ability to predict transitional flows was verified by 
comparing the results of the present procedure with the results of a well- 
established transitional boundary layer calculation procedure previously 
verified by an extensive comparison between theoretical predictions and 
experimental data (refs. 26, 28, 29, and 32). A comparison between steady- 
state transitional results of the present procedure, using the reduced set of 
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equations, and the results of the boundary layer procedure are presented in 
Fig. 8 . As can be seen, the results are nearly identical. The transition 
model utilized was virtually identical in both cases. 

Turbulent viscous flow calculations. - An additional comparison was made 
between the present procedure and the data of Karlsson (ref. 36) to assess the 
procedure's applicability to time -dependent turbulent flows. In ref. 36 
Karlsson investigated a boundary layer developing under an oscillating free- 
stream velocity given by 


u oe( x »t) “ u 0 ( x ) COSuJt 


(64) 


Karlsson represented the velocity field within the boundary layer by the 
series expansion 


u(x,y,t)«u(x,y)' 


u (l) cos^cosu>t - u (l) sin£smu>t + r 


(65) 


and measured u(x,y), cos <t> and u^sindi . In Eq. (65) r represents the 

higher order harmonics. Calculations corresponding to Karlsson 's data for the 
conditions 


u )/ 27 t - 1.0 cycles/sec 
u ( «/u 0 - 0.35 


(66) 


were carried out under two different sets of assumptions. The first calcula- 
tion assumed the boundary layer at time zero to be independent of the stream- 
wise coordinate, x; this is termed the similarity solution. The second 
calculation assumed that at an upstream station, the velocity satisfied an 
equilibrium turbulent boundary layer profile; i.e., 

u (y' ! ) ■ u t(y)(' ~ CGS<Jt ) ^ 67) 
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where u e (y) is the steady-state turbulent boundary layer profile. This time- 
dependent initial value solution was then allowed to develop as the flow 
proceeded dovnstream. The results of both calculations are presented in 
Fig. 9 where u(l)cos4> is represented by u^ and u^sini is represented by 
Uq''-*-). The agreement is moderately good. In both calculations the predicted 
mean velocity profile agrees well with data and the predicted first harmonics 
agree qualitatively with the data. Both calculations predict the rapid rise 
of the in-phase harmonic with distance from the wall and, although both 
calculations predict a velocity overshoot for the in-phase component, both 
disagree with Karlsson's measurements as to the magnitude and location of the 
overshoot. Similarly, agreement between theory and experiment for the out -of - 
phase component is good in the wall region but disagreement exists in the outer 
region. However, considering the present limits of the theoretical predic- 
tions, which either assume a similarity solution or a set of initial conditions 
at a given streamwise station in which no overshoot of the first in-phase 
harmonic is present and in which no first out -of -phase harmonic exists (see 
Eq. (67)), the agreement between theory and experiment seems satisfactory. 

Comparison between the full and reduced sets of equations. - As previously 
discussed, there are two options available within the viscous calculation 
procedure used in the present report, the option using the full stream 
function-vorticity relation, Eq. (4), and the option using the reduced stream 
function-vorticity relation, Eq. (19). The full relation represents an exact 
equation with no approximation tj the original Navier-Stokes equations, 
whereas the reduced equation requires the approximation 


±L 

dt 


<< 



( 68 ) 


where x is associated with the streamwise direction, y with the transverse 
direction, and u and v are velocity components in the x and y directions, 
respectively. The approximation obviously is valid in attached boundary layer 
type flows and, as shown by the calculations of Briley and McDonald (ref. 12), 
the approximation appears to be valid in relatively thin separation bubbles; 
however, at the initiation of the present study it was not obvious if the 
reduced set of equations leads to valid solutions in the stagnation region of 
the airfoil. Since the reduced set of equations are much more efficient to 
solve than the full set, they are used whenever possible and, therefore, a 
test comparison between the results of the full set and the reduced set was 
made for steady flow in the vicinity of an airfoil front stagnation point. 
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A comparison between the two results for a modified NACA 0012 airfoil at 
2.5 deg angle of attack is shown in Fig. 10. The Reynolds number for the 
calculation was 0.26 x 107 and a free -stream turbulence level of 0.01 was 
assumed. In Fig. 10 the normal distance and surface coordinate are normalized 
by the chord length. The normalized surface coordinate, s/c, is zero at the 
trailing edge and increases along the pressure surface to a value of 1.015 at 
the leading edge. The surface coordinate then continues to increase along the 
suction surface J o the trailing edge. The stagnation point is located well on 
the pressure side of the leading edge at s/c = 1.00. 

As shown in Fig. 10, the only significant difference in the results is 
reflected in the momentum thickness in the vicinity of the stagnation point. 
The full Navier-Stokes solution indicates a sharp increase in momentum 
thickness in the immediate area of the stagnation point due to the presence 
of a velocity overshoot, whereas the reduced Navier-Stokes equations solution 
gives a smooth variation of momentum thickness in the same region . The 
comparison indicates the validity of the solution from the reduced equations 
in the vicinity of the front stagnation point and, therefore, in subsequent 
viscous calculations the reduced set of equations was used exclusively. 


Predictions for Flow About Oscillating Airfoils 

General. - Weak -interact ion predictions were made for the flow field 
about a modified NACA 0012 airfoil for three different types of motion. The 
first calculation was for the flow about an airfoil oscillating in pitch 
about the quarter chord point. The flow conditions were based upon experi- 
mental data which showed that the airfoil did not encounter stall. The 
second calculation was made for the flow about an airfoil once again 
oscillating in pitch, however, in this second case the data showed the air- 
foil to be stalled over a large portion of the cycle. Finally, a third calcu- 
lation was made for an airfoil oscillating sinusoidally in heave; in this third 
case, the airfoil was stalled over a significant portion of the cycle. In all 
three cases the calculation was made by first predicting the inviscid flow 
field via the Giesing calculation procedure (ref. 10). This inviscid calcula- 
tion then serves as an outer edge boundary condition for the viscous flow 
calculation. The viscous flow itself is divided into subregions, as shown in 
Fig. 1. The viscous calculation has been described previously; in brief, the 
stagnation region is solved first and the results from the stagnation region 
serve as upstream boundary conditions for the next downstream region on both 
the suction and pressure surface. In this manner, the solution for each 
viscous region serves as a boundary condition for the next region. The mesn 
spacing in each region varied to allow adequate resolution of the local flow 
field, Typically, the streamwise mesh spacing is smallest in the leading 
edge separated region where Ax/c is of the order of 0 . 0025 . The largest 
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streamwise mesh spacing is used in the pressure side fully-turbulent region 
where Ax/c is roughly 0.0177. The distribution and number of grid points 
normal to the airfoil surface also vary from region to region so that optimum 
resolution of laminar boundary layers, transitional regions, and turbulent 
boundary layer sublayers is obtained. As an example of the total number of 
grid points used, the distribution grid points for each segment of the Case II 
airfoil are given below: 


Stagnation region 

Pressure side transition region 

Pressure side fully-turbulent region 

Suction side separated region 

Suction side fully-turbulent region 

Suction side fully-turbulent 
t-ailing edge region 

Total number of grid points 


30 x points 
34 x points 
34 x points 
30 x points 
34 x points 

34 x points 

5648 


by 24 y points 
by 26 y points 
by 33 y points 
by 26 y points 
by 30 y points 

by 23 y points 


Case I - unstalled airfoil oscillating in pitch. - The first time- 
dependent airfoil flew field calculation is for a modified NACA 0012 airfoil 
oscillating sinusoidally in pitch. The airfoil and motion chosen correspond 
to test point 3184.2 of ref. 3 for which the mean angle of attack is 7.76 deg, 
the amplitude of the sinusoidal oscillation is 5*24 deg, the dimensionless 
frequency, k, is 0.252, and the Reynolds number based upon chord length is 
0.26 x 107. The free -stream turbulence level is assumed to be one percent. 
Under static conditions, maximum lift for the modified airfoil at a chord 
Reynolds number of 0.26 x 10? occurs at 12.9 deg and the nonlinear portion of 
the lift -incidence curve begins at approximately 9 deg. Thus the airfoil is 
being investigated at conditions which slightly exceed the limit of static 
stall; however, from observation, the airfoil does not undergo any dynamic 
stall. Hie variation of the experimentally determined lift and moment 
coefficients compared with the theoretical predictions using the computer code 
of ref. 10 is presented in Fig. 11. The comparison between theory and 
experiment is co* 'idered good. 


A comparison between differential pressure coefficients, as predicted by 
the theory (ref. 10) and measured by experiment (ref. 3) , is presented in 
Figs. 12 through 14 and the surface pressure coefficient at several incidence 
angles is shown in Fig. 15 . As shown in Figs. 12 and 13 , the theoretical 
predictions are in good agreement with the data as might be expected for an 
unstalled airfoil. However, as shewn in Fig. 14, a large discrepancy does 
exist between theory and experiment at x/c » 0.91. The modified NACA 0012 
airfoil has a trailing edge tab (see Fig. l) and at the junction of this tab 
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and the airfoil skin there is a rapid change in surface curvature. This rapid 
change in curvature causes the potential flow program to seek a stagnation 
point and give a physically unrealistic pressure maximum a+ the tab location. 
For example, at an incidence of 7.5 deg the predicted pressure coefficients 
on the suction side trailing edge are as follows : 


x/c 

Cp 

0.67 

-0.18 

0.75 

-0.09 

0.82 

■HD . 06 

0.88 

+ 0.33 

0.93 

-KD.06 


The severe local maximum at the 88 percent chord location is physically 
unrealistic and is due to the sharp curvature in the airfoil surface at the 
juncture of the tab and the airfoil skin. In reality such a sharp maximum 
would not exist since the boundary layer displacement thickness would be 
expected to smooth the pressure distribution. Before the inviscid flow field 
calculated by the Giesing procedure is input to the viscous flow calculation, 
a three -point least -squares smoothing in time and a five -point least -squares 
smoothing in space is performed to insure a smooth variation of the outer edge 
viscous flow boundary conditions. However, additional smoothing was required 
in the region of the unrealistic pressure maximum. Prior to performing a 
least-squares smoothing of the inviscid velocity field, the local C p maximum 
(velocity minimum) is relieved somewhat through a two-point central averaging 
procedure. The averaging procedure consists of first obtaining the velocity 
at the 88 percent chord location as the average of the velocities just 
upstream and downstream of this location. The velocity at points upstream 
and downstream of the 88 percent chord point are then evaluated by marching 
both upstream and downstream away from the 88 percent location using the 
formula 


V M 


( 69 ) 


Considering j to be increasing in the direction of marching, it is noted that 
Vj_i is a previously averaged velocity. Equation (69) is used to smooth the 
velocity field between x/c * 0.80 and x/c = O.98; with this central averaging 
the only portion of the flew significantly modified is that in the immediate 
region of the 88 percent chord location. 
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The viscous flow field is initiated in the stagnation region using the 
reduced Navier-Stokes equations. As previously demonstrated, use of the 
reduced equations is quite justified on the basis of the comparisons with 
solutions to the complete Navier-Stokes equations and results in a consider- 
able reduction in computer run time. Use of the reduced equations still 
permits separation bubbles to appear as demanded by the physical constraints 
of the problem. The results predicted at both ends of this stagnation region 
are used as upstream conditions for two further segments, one on the suction 
surface and one on the pressure surface. The segment downstream of the 
stagnation region on the suction side of the airfoil contains the region where 
the leading edge separation bubble appears. This segment is then followed by 
the suction side fully-turbulent region which in itself may be divided into 
two or more segments. On the pressure side of the airfoil the stagnation 
region is followed by the segnent where boundary layer transition usually 
takes place. This region is then followed by the pressure side fully-turbulent 
region . 

Results of the calculations in the stagnation region are presented in 
Figs. 16 through 22. The calculation in this region is made assuming quasi - 
steady flow and the results are monitored to verify the validity of this 
approximation. Quasi -steady calculations are made in the stagnation region at 
discrete instants of time by using the instantaneous velocity distribution 
from the Giesing time -dependent inviscid flow procedure as a steady outer -edge 
boundary condition. The viscous calculation is carried out by assuming an 
initial viscous flow field and then letting the calculation march to a steady- 
state holding the outer edge velocity constant. The stagnation region calcu- 
lation is first made at a time ti in the cycle assuming the initial flew field 
to have a quadratic velocity profile. After this initial solution is 
converged, the next viscous calculation is made at a time t 2 in the cycle 
using the inviscid flow velocity distribution at time t£ and assuming as an 
initial guess that the viscous flow field at time t 2 is a scaled version of 
the viscous flow field at time t^. The scaled flow field gives the initial 
conditions for the viscous flow field and the calculation procedure is then 
allowed to march in time to a steady-state while holding the outer edge 
velocity constant. Thus the stagnate.. region is calculated as a series of 
quasi-steady solutions at selected points in the cycle by impulsively changing 
the outer edge inviscid flow velocity distribution and allowing the scaled 
viscous region to adjust to this new outer edge boundary condition. The time 
increment (tjj+i-tn) is typically 0.1 for the cycle. If the time required to 
adjust to this new condition is t ( the relaxation time), the ratio of 
T /(tn+l _ ta) * a typically of the order of 0.04. If d>i is the solution at t^ 
and $2 Is the solution at t 2 then t is typically taken to be the time 
required for (^-^) = 0.88 (dg'^i)* Ulus, since the solution approaches 
in a somewhat asymptotic manner, r can be regarded as two time constants of 
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an exponential decay from <$ ^ to <i> ^ or from the first steady solution to the 
second. For the purpose of these calculations the vorticity at each point in 
the flow field was monitored. The time constant, r, was defined as the time 
required for the vorticity at all points in the flow field to satisfy the 
criteria (<d-uji) = 0.88 (u^-u^) where uu i° the vorticity. Since the time 
required by the viscous layer to adjust to a new inviscid flow is so much 
smaller than the time required to change the inviscid flow, the viscous flow 
may be properly regarded as quasi -steady in the stagnation region. This 
quasi -steady nature of the viscous leading edge region has also been hypoth- 
esized by McCroskey (ref. 37) and Patay (ref. 38), based upon boundary layer 
calculations. 

The location of the stagnation point as a function of the Incidence angle 
is presented in Fig. 16. The surface coordinate, s, is zero at the trailing 
edge and increases as the coordinate travels along the pressure surface to the 
leading edge and then along the suction surface back to the trailing edge. 

The leading edge is located at s/c = 1.015 where c is the chord length. 
Although the viscous flow may be considered quasi-steady. Fig. 16 shows that 
the inviscid flow is certainly not quasi -steady since the front stagnation 
point location is not solely determined by the incidence angle but depends 
upon the flow time history. Predictions of skin friction coefficient in the 
stagnation region at selected incidence angles are presented in Fig. 17 where 
the distributions are presented such that time increases from the bottom of 
the figure to the top. Prediction of momentum thickness and displacement 
thickness are presented in Figs. 18 and 19, the location of the stagnation 
point is indicated by a circle. It should be noted that the displacement 
thickness, 6*, and the momentum thickness, 0, are defined by 

»*•{*('-£)<» < TO > 


® ‘ 4 * ^ ('■ ■&■) ** (?i) 


Since both edge velocity, u e , and the local velocity, u, are zero at the 
stagnation point, the Integrands are indeterminate; however, the integrals are 
not necessarily zero. It should be noted that a grid point was not located at 
the stagnation point in any calculation. As the flow proceeds away from the 
stagnation point in both the pressure and suction segments the integral 
thicknesses tend to increase due to the Influence of skin friction and tend to 
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decrease due to the influence of the very strong favorable pressure gradient. 
The results show the pressure gradient effect to be dominant on the suction 
side of the stagnation point and the skin friction effect to be dominant on 
the pressure side of the airfoil. It should be pointed out that the favorably 
pressure gradient on the suction side of the airfoil is greater than that on 
the pressure side and, therefore, the results obtained are not at all 
unreasonable. The variation of mixing length in the stagnation region is 
presented in Fig. 20. Near the stagnation point the nixing length grows due 
to the large amount of entrainment of free-streaur turbulence which acts as a 
source term on the turbulence kinetic energy equation, however, in this region 
the eddy viscosity is small compared to the molecular viscosity due to the 
relatively small transverse velocity gradients and the large viscous wall 
damping effect. Elsewhere in the stagnation region the mixing length is 
negligible. 

Predictions of velocity profiles in the stagnation region are presented 
in Figs. 21 and 22. Predictions of the profile at various incidence angles 
for a given location on the pressure side of the stagnation point are shewn 
in F' ’ . 21 and on the suction side of the stagnation point in Fig. 22 where 
Ua, is the free -stream approach velocity, ’’hen the incidence is at its maximum 
value, a = 13 deg, the stagnation point moves far down on the pressure side 
of the airfoil as demonstrated by the velocity plots. The profiles on the 
suction side of the airfoil are fuller than those on the pressure side 
reflecting the stronger favorable pressure gradients present on the suction 
surface. It should be noted that ali , -'ofiles appear to be laminar. 

The viscous flow in the pressure side segment containing the transition 
region is presented in Figs. 23 through 27. The variation of skin friction 
coefficient is presented in Fig. 23, where the skin friction coefficients at 
upstream locations in this region indicate laminar flow. As the flow 
progresses downstream it undergoes transition. As the incidence angle 
increases the transition is delayed due to the increased favorable pressure 
gradient and the movement of the stagnation point rearward on the pressure 
side. At the higher incidence angles transition is delayed' until near the end 
of this segment. Predictions of momentum and displacement thickness for the 
pressure side transition segment are presented in Figs. 24 and 25. The 
characteristic sharp increase in slope of the integrel quantities as a function 
of Ji 3 tance at the start of the transition zone is evident at the lower angles 
of attack. The variation of mixing length is presented in Fig. 26. The 
pressure side transition region also is calculated assuming a quasi -steady 
viscous flow field and in this region the ratio of viscous relaxation time to 
the inviscid flow change time is of the order of 0 . 15 , again indicating quasi- 
steady viscous flow. Rredictions of velocity profiles at a constant stream- 
wise station in the pressure side transition region are shewn in Fig. 27. The 
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streamwise station was chosen to be near the end of the segment and the 
velocity profiles show the fully-turbulent, transitional, and laminar 
character of the flow as the transition zone moves rearward and then forward 
as the incidence angle increases and decreases, respectively. Hie results on 
the fully-turbulent pressure segment are presented in Figs. 28 through 32 
which again show skin friction coefficient, displacement thickness, momentum 
thickness, mixing length distributions, and velocity profiles at specified 
angles of attack. However, in this segment the quasi-steady approximation is 
not valid and the flow in this region was predicted as time -dependent and 
viscous. 

The predictions of the flow field for the segment containing the .parated 
region on the suction side of the airfoil are presented in Figs. 33 nrough 
38. At lower incidence angles no bubble is present, however, as the incidence 
angle increases a bubble does appear. The position of the bubble and the 
resulting streamline pattern at seven specified incidence angles are presented 
in Fig. 33. No bubble is present at a = 2.5 deg. The theory predicts a 
bubble to appear in the vicinity of the leading edge f the airfoil by the 
time the incidence is equal to 6.11 deg. As the incidence angle increases the 
bubble increases in size and moves upstream. A further increase in incidence 
leads to a -further upstream movement and contraction of the bubble. As the 
incidence angle is decreased the bubble moves backward on the airfoil, 
increases in size, then decreases in size, and finally disappears. The bubble 
behavior predicted by the theory corresponds to experimental observation for 
hover tests of Velkoff, Blaser, and Jones (ref. 39) and was hypothesized by 
Ham (ref. k) based upon various steady-state data and the time -dependent 
data of Isogai (ref. 40). Ham attributed the dealy in stall 
under dynamic conditions to the retardation of the bubble reattachment point 
with incidence angle under dynamic conditions. In the present study the 
separated viscous flow region was calculated under the viscous quasi -steady 
assumptions. The ratio of viscous relaxation time to inviscid flow change 
time was of the order of 0.15, thus indicating that the quasi -steady viscous 
flow assumption is valid. Therefore, the present results indicate that, if 
Ham's hypothesis is to be believed, the delay in the movement of the 
reattachment point is due to the time -dependent nature of the inviscid flow. 
Distributions of skin friction, momentum thickness, and displacement thickness 
for the segment containing the separated region are presented in Figs. 3*+ 
through 36. The skin friction plot indicates clearly the length of the 
separated region. Hie streamwise variation of mixing length for this segment 
is presented in Fig. 37. The turbulence kinetic energy model has been used 
to predict transition and the turbulence field in this region. As shown in 
Fig. 37, the predicted bubbles are transitional. At any given instant the 
bubble behavior predicted by the present theory and presented in Figs. 33 
through J7 is in qualitative agreement with McCullough and Gault's (ref. l) 
description of steady-state leading edge bubbles. 
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The velocity profiles in the suction side separated region at selected 
instantaneous incidence angles are shown in Fig. 38 at a fixed streamwise 
location. Moving from left to right, the initial profile, a = 2.5 deg, shows 
that the flow is approaching separation, however, at the low incidence no 
leading edge separation occurs. The second profile, a= 6.11 deg, shows the 
initiation of flow separation as the bubble appears on the airfoil (see Fig. 
33). The third profile, a = 9.35 deg, clearly reflects the forward motion of 
the bubble as the profile shows a larger reversed flow region than the previous 
profile. As the angle of attack continues to increase, the separation bubble 
becomes smaller and moves even further forward cm the airfoil. The fourth 
profile shown in Fig. 38 at a = 13 deg shows that at this point in the cycle 
the entire separation bubble has moved upstream of s/c = 1.057. The shape of 
the profile indicates the transitional behavior of the boundary layer in the 
area of flow reattachment. The last three profiles shewn in Fig. 38 then show 
the rearward movement, growth, and final disappearance of the separation 
bubble as the incidence returns to its minimum value. 

An examination of the suction side inviscid velocity distribution 
predicted by the Giesing calculation procedure shows a strong adverse pressure 
distribution followed by a strong favorable pressure distribution ir. the 
region of the tab which is located at the trailing edge of the airfoil (see 
Fig. 1). Even after the smoothing procedure of Eq. ( 69 ) was carried out it 
was felt that adequate resolution in the vicinity of the tab would require 
more grid points than could be accommodated in core storage if the entire 
suction fully -turbulent region were done simultaneously and, therefore, the 
suction side fully-turbulent regicn itself was divided into two regions both 
of which were calculated as unsteady regions. The first region is between the 
10 percent chord station and the 80 percent chord station and the second 
region is between the 80 percent chord station and the trailing edge. 
Calculations in the first of these regions are shown at specified angles of 
attack in Figs. 39 through 42. The calculations on the aft section of the 
airfoil are presented in Figs. 43 through 46. Although the suction side 
trailing edge region appears to approach separation (defined as the appearance 
of a region of reversed flow), separation is never reached. As shewn in 
Figs. 42 and 46, the flew remains turbulent in these regions. Velocity 
profiles in both fully-turbulent regions are shewn in Fig. 47. 

Case II - the stalled airfoil oscillating in pitch. - The second set of 
calculations made are again for an airfoil oscillating in pitch. This case 
corresponds to test point 3171.4 of ref. 3 for which the mean incidence is 
12 . 5,3 deg, the amplitude of oscillation is 5.39 deg, the reduced frequency, k, 
is 0.242, and the Reynolds number based upon chord Is 0.47 x 10?. Under static 
conditions, the maximum lift for a Reynolds number of 0.47 x 10? occurs at 
9.5 deg, therefore, the airfoil is operating well above the static stall limit 
over most of its cycle. 
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A comparison between theoretica 1 v predicted and experimentally measured 
lift and moment coefficients is prest *d in Fig. 48. As can be seen, serious 
disagreement exists between theory and experiment, particularly for the moment 
coefficient at high incidence. Throughout the entire cycle the predicted lift 
coefficient is considerably higher than that measured. This comparison 
indicates that in contrast to Case I which is not in stall, the Case II air- 
foil is in stall. The differential pressure coefficients, shown in Figs. 49 
through 51 » also indicate the airfoil is experiencing a more extensive stall 
at the higher angles of attack. Surface pressure coefficients are shown in 
Fig. 52. A comparison between Figs. 50 and 51 for theoretical predictions at 
x/c = 0.88 and 0.91 indicates the very poor local prediction at the tab. As 
in the previous calculations, the velocity distribution in the region of the 
locally predicted pressure maximum was smoothed using the averaging procedure 
outlined in Eq. ( 69 ). This distribution was then least -squared and input 
into the viscous flow calculation procedure. 

The viscous calculations for the Oase II airfoil divide the viscous flow 
into six regions, including a stagnation region, a pressure side transition 
region, a pressure side fully-turbulent region, and a suction side separation 
region, as shown in Fig. 1. In addition, the airfoil suction side fully- 
turbulent region is divided into two subregions. As for the Case I airfoil 
the coordinate, s, represents the distance along the surface. Results for the 
stagnation region are presented in Figs. 53 through 57- As in Case I, it was 
found that the stagnation region, the pressure side transition region, and the 
suction side separated region are quasi -steady. In these regions the time 
required for a viscous solution tc adjust from a set of outer boundary 
conditions at time, tj_, to a new set of outer boundary conditions at a time, 
t^ + At, is much less than At. The fact that the calculated viscous 
response time scale is much smaller than the outer flew inviscid time scale 
demonstrates that the flow in these regions is quasi -steady. 

The location of the stagnation point as a function of the incidence angle, 
shown in Fig. 53, indicates that the stagnation point moves between s/c = 

O .926 and s/c = O. 983 . Since the airfoil nose is at s/c = 1.015, the stagna- 
tion point is always on the underside of the airfoil. The variations of skin 
friction in the stagnation region at various instantaneour incidence angles are 
shown in Fig. 5 1 '-. In this region the skin friction coefficient shows its 
expected high value in the vicinity of the stagnation point and then rapidly 
decreases. The variation in momentum thickness and displacement thickness at 
various instantaneous angles of attack is presented in Figs. 55 and 56 . The 
location of the stagnation point is denoted in Figs. 55 and 56 by an open 
circle. The predicted mixing length distribution, jt a /6, at several incidence 
angles is presented in Fig. 57. In the region of the stagnation point the 
mixing length is fairly large due to a large amount of entrainment of free 
stream turbulence. However, in the imnediate vicinity of the stagnation point 
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the transverse velocity gradients are small and the viscous damping is large 
leading to the eddy viscosity being much smaller than the laminar viscosity. 
Thus the flow acts laminar. As the flow proceeds away from the stagnation 
point, the mixing length decreases but eventually may increase depending on the 
turbulent energy balance. 

Predictions of the flow in the pressure side transition region are 
presented in Figs. 58 through 6l. The skin friction distributions of Fig. 58 
and the mixing length distributions of Fig. 6l show the flow to be undergoing 
transition from the laminar to turbulent state. The transitional nature of 
the viscous flow is also indicated by the change in slope of the momentum 
thickness and displacement thickness plots of Figs. 59 and 60. The solutions 
for the fully-turbulent pressure side region are presented in Figs. 62 through 
65. In this region the flew is unsteady and no quasi-steady assumption is 
used. The flow remains fully-turbulent here as demonstrated by the mixing 
length distributions presented in Fig. 65. 

The segment downstream of the stagnation segment on the suction side of 
the airfoil is termed the separation region. In this small region extending 
over 6 percent chord, there is a large adverse pressure gradient and *+ high 
incidence angles a leading edge separation bubble is likely to exist. The 
viscous flow predictions in the separated region are presented in Figs. 66 
through 70. A plot of the streamlines at various instantaneous incidence 
angles is presented in Fig. 66. The results show that a leading edge 
separation bubble is present through all of the cycle. The cycle proceeds 
from 0= 7.17, 11.92, 15.43, 17.92, 15.27, 11.74, and 7.17 in Fig. 66. 

Although not shown, there does exist a small portion of the cycle for a larger 
than its minimum value when no bubble is present. As shewn in Fig. 66, the 
bubble demonstrates the expected behavior of appearing as a relatively large 
bubble, moving forward as incidence angle increases and moving aft as incidence 
angle decreases. There is also a tendency for the bubble to shorten at the 
highest incidence angles. 

As previously discussed, this predicted behavior is confirmed by 
experimental evidence (refs 1, 38 through 40). Predictions of skin friction 
coefficients are presented in Fig. 67 and of momentum and displacement thick- 
ness in Figs. 68 and 69. The predicted mixing length distributions of Fig. 70 
indicate that the bubbles are transitional as the predicted mixing length is 
beginning to increase at stations at which the bubble exists. It should be 
noted that, although the experimental evidence indicates the airfoil is in 
stall through much of its cycle (Figs. 48 through 50), the predicted leading 
edge bubble behavior is not significantly different from that predicted in the 
unstalled airfoil, Case I. The major difference in the leading edge bubble 
region is that for the Case I airfoil the bubble is present over only approxi- 
mately one -half the cycle, whereas for the present Case II airfoil the bubble 
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is present over nearly all of the cycle. However, as will be shown 
subsequently, in the second suction side fully-turbulent region significant 
differences between the Case I and Case II calculations appear. 

The suction side fully-turbulent region is divided into two segments, 
the first extending from the 5 percent chord station to the 30 percent chord 
station, and the second extending from the 30 percent chord station to the 
trailing edge. The predicted results in the first suction side fully- 
turbulent region are presented in Figs. 71 through 74. As can be seen, the 
flow appears to be well-behaved in this segment. 

Although the first segment of the suction side Case II fully-turbulent 
region gives predictions similar to the predictions of the Case I airfoil, the 
second segment shews a very much different behavior between the Case II and 
Case I airfoils. It should be recalled that the Case I airfoil was determined 
experimentally not to ' e in stall, whereas the experimental data shows the 
Case II airfoil to be in stall through a significant portion of its cycle. As 
shown in Figs. 1+3 through 46, which present the suction side trailing-edge 
region for the Case I airfoil, the viscous flow approaches separation in the 
vicinity of the trailing edge but trailing edge separation does not occur. It 
may be expected that for an airfoil at higher incidence angle, such as the 
Case II airfoil, separation would occur in the trailing edge region. In order 
to understand the results presented for the Case II trailing edge region, it 
is helpful to review how the calculation in this region is made. 

The trailing edge region is a region in which time -dependent effects are 
important and, therefore, a time -dependent solution of the reduced Navier- 
Stokes equations is obtained. Since the segment is solved as a time -dependent 
flow field, it is necessary to specify an initial flew field at time to- This 
initial flew field is set by assuming that at the upstream boundary of the 
segment being calculated the flow variables are identical to those at the 
downstream boundary of the previous upstream segment at the same time t Q . The 
flow field is then assumed similar in the streamwise coordinate so that at any 
streamwise station at the initial time, t 0 , the flow variables are scaled 
distributions of the upstream flew variables. The distributions are scaled 
so as to match the outer edge velocity imposed by the time -dependent inviscid 
flow calculation procedure. Hie flow field in the segment being investigated 
is then calculated by letting the assumed initial flow field develop in time 
according to the governing equations. At each time, t n , the upstream boundary 
conditions are set equal to the conditions calculated at the downstream 
boundary of the previous segment at time, t n ; the outer edge boundary condi- 
tions correspond to the conditions calculated at time t n by the inviscid flow 
calculation procedure. Therefore, in summary, an initial flew field is 
assumed, time -dependent upstream and outer edge boundary conditions obtained 
from other calculations are imposed and the governing equations are solved. 


43 


Since the outer flow is cyclic, the effect of the initial conditions in 
general should disappear and a cyclic viscous solution ought to be calculated. 
For previous calculations the initial conditions in the viscous region 
disappear within less than one-quarter of a cycle; that is, the viscous flow 
solution became cyclic within one and one -quarter cycles of time integration. 
This was confirmed by carrying out calculations of the time -dependent region 
over approximately one and one -half cycles. 

Calculations for the second segment of the Case II fully-turbulent region 
are presented in Figs. 75 through 77. In contrast to the Case I calculations 
in which a cyclic nonseparated viscous region was calculated, the Case II 
calculations show a large separation bubble to appear. This behavior is 
demonstrated in Fig. 75 which shows skin friction coefficient at various 
incidence angles. The coefficient for the initial profile is shown at the 
bottom of Fig. 75 for an incidence of 7.17 deg. As shown in the remaining 
plots of Fig. 75, the skin friction in the trailing edge region drops rapidly 
as time and angle of incidence increase until at a = 17-92 deg a significant 
separation region (indicated by negative skin friction) is present. The 
separated region continues to increase in extent even after the incidence 
reaches a maximum and decreases. For example, at 15.27 deg the separated 
region extends over approximately 50 percent of the suction side of the air- 
foil. The calculation was continued back to the minimum angle of attack and 
no tendency for the bubble to disappear was noted. Prediction of momentum 
thickness and displacement thickness are presented in Figs. 76 and 77, and 
predicted streamlines are shown in Fig. 78. An examination of the results of 
Figs. 75 through 77 clearly indicates that the bubble is being prevented from 
moving upstream only by the upstream initial condition being imposed upon the 
segment. Hie dimensionless mixing length originally remained at about 0.09 
but as time increased grew to a value of approximately 0.40 in the separated 
flow regime. In addition, the bubble has, of course, grown to such an extent 
that the weak-interaction assumption of the viscous layer not affecting the 
pressure distribution is being violated. The effect of the bubble growth would 
relax the pressure gradient and in a strong -interaction solution which 
includes the mutual interaction between the viscous and inviscid flow fields 
an equilibrium may be reached and the bubble may stop growing. Within the 
limits of weak-interaction theory the indication is that the leading edge 
bubbles in the Case I unstalled and Case II stalled airfoils behave very 
similarly; no dramatic flow phenomena in the leading edge region are predicted 
in the stalled case. However, the trailing edge regions behave quite 
differently. In the unstalled case, Case I, the flow remains attached; in the 
stalled case, Case II, a trailing edge bubble is formed which moves rapidly 
upstream and separates the flew over a significant portion of the airfoil 
suction surface. This behavior suggests a possible stall mechanism for this 
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type of airfoil in which a trailing edge bubble is formed and then rapidly 
moves toward the leading edge as incidence is increased until a large portion 
of the airfoil is separated and stall occurs. 

Case III - the stalled airfoil oscillating in heave. - The final set of 
oscillating airfoil calculations are for an airfoil oscillating sinusoidally 
in heave. The airfoil and motion originally chosen correspond to test point 
3090.2 of ref. 3 for which the incidence angle in absence of vertical motion 
is equal to 12.36 deg, the magnitude of the vertical oscillation is O.306 
based on the semichord, the reduced frequency, k, is 0.242, and the Reynolds 
number based upon chord is 0.26 x 10?. Under static conditions, the maximum 
lift for a Reynolds number of 0.26 x 107 occurs at 12.9 deg and the nonlinear 
portion of the lift -incidence curve begins at approximately 9 deg. Thus the 
airfoil is being investigated at conditions which exceed the limits of static 
stall . 

Theoretical prediction of the lift and moment coefficients made using 
the theory of Giesing (ref. 10) are compared with the experimentally measured 
values of ref. 3, in Fig. 79* The predictions are indicated by the chain line 
and termed Run A. As can be seen, there is considerable discrepancy between 
theory and experiment, particularly for the moment coefficient at large 
negative values of normalized transverse location. The predicted lift 
coefficient is higher than the measured value throughout the entire cycle. 
However, some discrepancies do appear between the data of Liiva (see ref. 4l) 
and that of Halfman, Johnson, and Haley (ref. 42) and Rainey (ref. 43). In all 
three cited references experimental investigations were made for airfoils 
oscillating in heave. The data of both refs. 42 and 43 showed that the airfoil 
dissipates work to the surrounding air and, therefore, the airfoil is stable in 
the bending mode. The data of Liiva (refs. 3 and 4l) on the other hand 
disagrees and shows regions of bending mode instability to exist under the 
influence of heave (ref. 4l). Although this stability discrepancy is a 
discrepancy in the aeroelastic response, it may result from a discrepancy in 
the flow fields present about the airfoils. In view of these discrepancies 
between the Liiva data (ref. 4l) and that of Halfman (ref. 42) and Rainey 
(ref. 43) and the fact that the airfoil under the stated conditions only 
appeared to be stalled at the high negative values of the transverse coordi- 
nate, a second inviscid calculation, Run B, was made in which all motion 
parameters remained the same except the vertical oscillation was set at 0.612 
based on semichord. In this second case, Run B, the maximum instantaneous 
incidence angle is approximately 2C deg as opposed to a maximum of 16 deg in 
the original case and thus the second case, Run B, is a case in which the air- 
foil is much more likely to be in the region of stall over a significant 
portion of the motion cycle. Since it was desired to calculate the viscous 
flow field under conditions for which stall does occur and since the heave 
data of refs. 3 and 4l showed discrepancies with that of refs. 42 and 43, 
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the second inviscid run was used for the viscous calculations to insure 
calculation of a stalled airfoil. Theoretical predictions of lift and moment 
coefficient made using the Giesing procedure (ref. 10 ) are shown by * ">* solid 
line for Run B in Fig. 79* The differential pressure coefficients are shown 
in Figs. 80 through 82. Surface pressure coefficients for Run B are presented 
in Fig. 83. A comparison between Figs. 8l and 82 for theoretical predictions 
at x/c = 0.88 and 0.91 indicate the very poor local prediction at the airfoil 
tab which was also present for both the Case I and Case II airfoils. As in 
the previous calculations, the velocity distribution in the region of the 
locally predicted pressure maximum was smoothed using the averaging procedure 
of Eq. ( 69 ). This distribution was then least-squares smoothed and input into 
the viscous flow calculation procedure. 

The viscous calculations for the Case III airfoil divide the viscous flow 
field into only five regions instead of the six regions used for the Case I 
and Case II airfoils. The division is as follows: a stagnation region, a 

pressure side transition region, a pressure side fully-turbulent region, a 
suction side separation region, and a suction side fully-turbulent region. In 
the calculations made for the Case I and Case II airfoils, the suction side 
fully-turbulent region was further divided into two segments. However, due to 
the results of the calculations for the Case II airfoil in which a large 
trailing edge separation bubble was prevented from moving forward due to the 
location of the segment boundary, the suction side fully-turbulent region was 
treated as a single segment for the Case III airfoil. As in Case I and Case 
II, it was again found that the stagnation region, the pressure side transi- 
tion region, and the suction side separated region were quasi -steady. Thus 
the flow in these regions was calculated using the quasi -steady viscous flow 
assumption although the imposed outer edge velocity distribution was that 
obtained from the time -dependent inviscid flow computation. 

The results of the calculation for the stagnation region are presented in 
Figs. 84 through 88. The variation of skin friction in the stagnation region 
at various instantaneous transverse positions is shown in Fig. 85 . In this 
region the skin friction coefficient shows its expected high value in the 
vicinity of the stagnation point and then rapidly decreases. The variation in 
momentum 'thickness and displacement thickness at various instantaneous 
transverse positions is presented in Figs. 86 and 87 . The stagnation point is 
indicated by an open circle. The predicted mixing length distribution, l 0 Jb, 
at several transverse positions is presented in Fig. 88. 

Predictions of the flow in the pressure side transition region are 
presented in Figs. 89 through 92. The skin friction distributions of Fig. 89 
and the mixing length distributions of Fig. 92 show the flow to be undergoing 
transition at ail but two transverse positions. The transitional nature of 
the viscous flow field is also indicated by the change in slope of the 
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momentum thickness and displacement thickness plcts of Figs. 90 and 91* The 
solutions for the fully-turbulent pressure side region are presented in Figs. 

93 through 96. In this region the flow is unsteady and no quasi-steady 
assumption is used. Hie flow remains fully-turbulent here, as demonstrated by 
the mixing length distributions presented in Fig. 96. 

Hie Case III separation region extends over 7 percent chord where there 
is a large adverse pressure gradient and due to the high incidence angle 
(recall that the geometric incidence angle is constant for this case) a lead- 
ing edge separation bubble is likely to exist. Hie viscous flow predictions 
in the separated region are presented in Figs. 97 through 101. A plot of the 
streamlines at various instantaneous transverse locations is presented in 
Fig. 97* Hie results shew that a leading edge separation bubble is present 
through the complete cycle. Hie cycle proceeds from Y/Ymnv = -1.0, -0.309, 
+0.309, +1.0, -K).309> -0.309, and -1.0 in Fig. 97 with time increasing from 
bottom to top. Note that initially the airfoil is in its lowest position and 
the separation bubble has moved forward on the airfoil. Then, as the airfoil 
moves upward, the bubble moves back and grows in size. As the airfoil begins 
to move downward, the separation bubble becomes shorter and moves forward 
again. This behavior is explained by considering the relative incidence angle 
of the airfoil. Although the incidence angle (measured with respect to a 
constant reference direction) is constant, as the airfoil moves upward its 
velocity is added in a vector sense to the velocity of the free stream and the 
resultant velocity yields a lower effective incidence angle. Likewise, when 
the airfoil moves downward, the effective incidence angle is increased. With 
this behavior in mind, the predicted results appear to be in agreement with 
what is observed experimentally for airfoils oscillating in pitch (refs. 1, 
and 38 through 40). Predictions of skin friction coefficients are presented 
in Fig. 98 and of momentum and displacement thickness in Figs. 99 and 100. 

The predicted mixing length distributions of Fig. 101 indicate that the 
bubbles are transitional as the predicted mixing length is beginning to 
increase at station^ in which the bubble exists. 

Hie suction side fully-turbulent region of the Case III airfoil was 
treated as one segment rather than the two used for the Case T and Case II 
airfoils. This was accomplished by using variable mesh spacing in the stream- 
wise direction. The streamwise mesh spacing varied so that a relatively 
tight mesh was used at the upstream boundary of the segment. The mesh spacing 
then increased until the region of the tab was approached. Near the tab the 
spacing was again decreased. Hie calculation was performed in this manner so 
that any separation bubble which appeared would have a segment of nearly 
90 percent of the airfoil chord over which it could grew without encountering 
the upstream boundary of the segment as it did in the Case II calculation. 

Hie results for this segment are shown in Figs. 102 through 104. As in the Case 
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II calculation, a cold start initial condition is shown at the bottom of Figs. 
102 through 104. As discussed previously, this cold start initial condition 
is obtained by specifying the stream function and vorticity profiles at the 
upstream boundary and then scaling these profiles at each ctreamwise station 
to match the edge velocity. As the calculation progresses in time a large 
separation bubble appears almost immediately and continues to grow. This is 
shown in the plots of skin friction (Fig. 102) where, for example, when the 
airfoil is at a transverse position of +1.00 the separated region is almost 
50 percent of the chord in length. When the airfoil begins to move downward 
and the effective incidence angle becomes greater, the separated region grows 
rapidly until at Y/Y rr^v = +0.8 it covers nearly 70 percent of the airfoil, as 
indicated by the top plot of skin friction coefficient shown in Fig. 102, at 
which time the calculation was terminated. The momentum thickness and dis- 
placement thickness are shown in Figs. 103 and 104. Figure 103 shows that the 
displacement thickness at the termination of the run is extremely large and 
the assumption of weak-interaction is severely violated. For this reason and 
due to the similar behavior of the Case II calculation, no attempt was made to 
carry out the computation for the remainder of the cycle, although the 
calculation was still numerically stable. 


DISCUSSION AND CONCLUSIONS 


A weak-interaction solution for the problem of the flow field about an 
airfoil in arbitrary unsteady motion has been developed by combining an 
unsteady nonlinear potential flow computer code (ref. 10) with a finite - 
difference viscous flow computer code (ref. 12). The potential flow proce- 
dure serves to predict an inviscid flow field about the airfoil and this 
inviscid flow field is input into the viscous procedure as an outer edge 
boundary condition for the viscous layer. The viscous development is then 
predicted under the influence of the applied inviscid flow field using the 
weak-interaction assumption that the viscous flow does not significantly 
influence the outer inviscid flow field. The weak-interaction assumption is 
valid as long as the viscous displacement thickness remains small compared to 
the airfoil thickness. However, when the displacement thickness becomes large 
and significantly modifies the inviscid pressure distribution, such as in a 
region of significant boundary layer separation, the weak-interaction theory 
is invalid and accurate predictions of the flow field under these conditions 
requires a theory which recognizes the mutual interaction between the viscous 
inner and nominally inviscid outer flow fields. Such a strong-ir.teraction 
calculation procedure could be developed by an extension of a successful 
weak-interaction procedure in which an inner viscous solution such as the 
viscous solution of the present report is coupled to an inviscid outer solu- 
tion. The coupling would require continuity of flow angle along the line 
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joining these solutions. Alternatively, the entire flow field could be 
solved by the Navier -Stokes equations thus avoiding the problem of coupling 
twc different solutions in two regions of the flow. 

Although a weak-interaction solution is limited in applicability to flow 
situations in which the viscous displacement thickness does not significantly 
affect the inviscid pressure distribution, weak-interaction solutions should 
accurately predict airfoil flow fields if no significant regions of separation 
are present. In this regard the weak-interaction solution can give a 
quantitative picture of viscous flow phenomena such as demonstrated by the 
separation bubble calculation presented previously. In addition, the procedure 
should be able to predict incipient stall since when incipient stall occurs 
the separated region is still small enough to allow the weak-interaction 
assumption to be valid. It should be noted that in their study of transitional 
separation bubbles Briley and McDonald (ref. 12) included a strong-interaction 
viscous correction to the inviscid pressure field. However, this was a 
linearized correction and is only valid for thin separated regions. Therefore, 
it could not be validly applied to the thick trailing edge separated regions 
computed in the present effort. 

Three airfoil calculations have been presented; the first of these, 

Case I, corresponds to an airfoil experimentally found not to exhibit any 
characteristics of stall, and the second and third calculations, Case II and 
III, were for airfoils which are in stall over a significant portion of the 
motion cycle. The viscous calculations in the first unstalled case differed 
considerably from the calculations in the second and third stalled cases. In 
the unstalled case a well-behaved viscous flew was found to exist over the 
entire airfoil, whereas in the stalled cases significant separated regions 
appeared along the trailing section of the airfoil. Leading edge separation 
appears on the suction surface in all cases; however, in both the stalled and 
unstalled calculations the flow in the leading edge separation region soon 
undergoes transition, reattaches, and forms a well-behaved bubble. Thus, no 
qualitative differences in the leading edge separation bubbles are predicted 
between the experimentally observed unstalled and stalled airfoils. 

The predicted leading edge bubble behavior indicates a possible mechanism 
of leading edge stall. As the airfoil increases in incidence, a bubble 
appears in the leading edge region. In each case the boundary layer is laminar 
at the leading edge bubble separation point and undergoes transition to the 
turbulent state. The calculations show an increase in incidence to be 
accompanied by a forward movement and a shortening of the bubble, as has been 
deduced from experimental data by Velkhoff, Blaser, and Jones (ref. 39) • 

Isogai (ref. 4o), and McCullough and Gault (ref. l). In each case the forward 


49 



movement of the separation point with increasing incidence angle is 
accompanied by a forward movement in transition location and, hence, a 
subsequent forward movement in reattachment. The net result being a predicted 
shortening of the separation bubble as incidence increases. 

Whether or not the bubble undergoes transition is determined by both the 
applied pressure gradient and the free-stream disturbance level. If the 
bubble does not undergo transition, it is expected that the bubble would not 
reattach and leading edge stall would occur. Within the limits of weak- 
interaction theory this would correspond to a viscous solution in which the 
leading -edge bubble were to grow very rapidly leading to a large separated 
region. In the present calculations transition always occured soon after 
separation; however, transition is a strong function of the free-stream 
turbulence level (ref. 26) and, if a low enough turbulence level were assumed, 
transition would be expected to be delayed and the leading edge bubble could, 
depending on local conditions, gr <4 rapidly leading to eventual leading edge 
stall. 

The major difference in the viscous calculations betwen the stalled and 
unstalled airfoils occurs in the trailing edge suction side region. In Case I, 
the calculation for the air foil which has been determined experimentally not 
to be in stall, the suction side trailing edge segment is a well-behaved 
viscous flow region. The viscous layer approaches separation at the junction 
of the airfoil skin and the trailing edge tab but separation does not occur. 

In Case II and Case III, both of which have been determined experimentally to 
be in stall, a large separated region does appear along the suction side 
trailing segment. Due to the absence of a strong -interaction mechanism for 
alleviating *e pressure distribution, the separated region grows uncontrolled. 
In Case II the trailing edge separated region encompasses over 50 percent of 
the airfoil suction surface and is prevented from becoming larger only by the 
location of the segment boundary. In Case III the trailing edge separated 
region grew to 70 percent of the airfoil surface and the displacement thick- 
ness reached approximately 25 percent of the airfoil chord when the calcula- 
tion was terminated. The appearance of this large separated region is 
interpreted as indicative of stall. Thus the viscous behavior in the 
unstalled and stalled cases is significantly different. In the unstalled 
case, Case I, the flow remains attached in the trailing edge region, whereas 
in the stalled cases, Case II and Case III, a trailing edge bubble is formed 
which moves rapidly upstream and separates the flow over a significant portion 
of the airfoil suction surface. Uais behavior suggests a possible stall 
mechanism for this type of airfoil and motion in which a trailing edge bubble 
is formed and then rapidly moves toward the leading edge until a large portion 
of the airfoil is separated and stall occurs. This mechanism could cause 
trailing edge stall or could modify the overall pressure distribution about 
the airfoil in such a manner as to cause the leading edge bubble to fail to 
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reattach; this failure of the leading edge bubble to reattach would cause the 
airfoil to go into leading edge stall. In regard to the problem of whether 
the NACA 0012 airfoil actually undergoes leading edge or trailing edge stall, 
it has been pointed out by Ericsson and Reding (ref. 6) that under steady 
conditions the NACA 0012 airfoil is a prime candidate for transition between 
stall types with variation in chord Reynolds number; according to ref. 6, the 
stall shifts from the leading edge type to the trailing edge type t a Reynolds 
number based upon chord of approximately 0.6 x 107. Therefore, th* NACA 0012 
airfoil could be expected to exhibit either leading edge or trailing edge stall 
or ever, a combination of the two. 

Hie viscous flow calculation procedure divides the entire viscous region 
into several subregions and cm examination of the results in each of these 
subregions indicates the degree of sophistication required ir. each portion of 
the viscous flow field. The calculations indicate the leading edge region to 
be quasi-steady; i.e., the viscous shear layer adjusts to changes to the outer 
edge boundary conditions in a time scale much shorter than the time s-.ale of 
the outer flow. This quasi-steady conclusion has also been reached by 
McCroskey (ref. 37 ) and Patay (ref. 38 ) based upon boundary layer calculations. 
When the leading edge region is quasi -steady and when no viscous flow 
separation appears, the viscous region can (and should) be calculated through 
a finite -difference steady-state boundary layer procedure. In general, 
finite -difference, steady-state boundary layer procedures would be expected 
to be considerably more rapid in terms of computer time than a solution of the 
full Navier-Stokes equations or than the solution of the 'reduced' Navier- 
Stokes equations used in the present effort. It should be pointed out that 
the reduced Navier-Stokes equations are equivalent to a set of time “dependent 
boundary layer equations with the addition of a streamwise diffusion term and 
solve the steady-state problem by assuming an initial flew field and then 
allowing the flow field to develop in time under steady-state boundary 
conditions. However, the present procedure for solving the reduced equations 
is highly competitive with the more usual time -dependent boundary layer codes. 

The pressure side of ^he airfoil downstream of the stagnation region 
consists of a transition region and a fully -turbulent region. The transition 
region has l fen found to be quasi -steady for the airfoils investigated under 
the pie sent effort. Therefore, as for the stagnation region, a steady-state 
boundary layer procedure may be preferable to the asymptotic time solution of 
the reduced set of equations used in the present study. However , the best 
available transition model must be enibodled in any computation procedure for 
the pressure side transition region. The fully-turbulent pressure side region, 
was found to be unsteady and it is felt that the procedure which was used 
(i.e., the solution of the reduced set of equations) was competitive with any 
presently available alternative procedure. 
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On the suction surface the stagnation region is followed by the region 
where a leading edge separation bubble is expected to appear. This region 
contains a very complex flow field which exhibits bolt eparation and 
transition. Any procedure used in this region should be capable of accurately 
predicting both phenomena. As has been shewn by the predated separation 
bubble behavior, the present procedure appears to predict the experimentally 
observed physical behavior of leading edge separation bubbles. Based upon the 
excellent qualitative predictions of the present procedure and the extreme 
complexity of the flow field in this region, it does not appear that any 
simplified analysis should be used to predict the leading edge separation 
bubble behavior. The leading edge separation region is followed by a suevien 
side unsteady, fully-turbulent flow region. Once again it is felt that the 
procedure for solving the reduced set of equations, which was used in this 
region, should be competitive with any pre ently available alternate procedure. 

In summary, the major modifications to the current procedure which would 
be recommended in a weak -interaction solution are in the stagnation region and 
the pressure side transition region. In both regions a savings in computer 
running time could be gained by replacing the present solution procedure with 
a steady-state finite -difference boundary layer procedure using the best 
available transition model. In all other regions the procedures used in the 
present effort appear to contain a necessary and sufficient amount of 
sophistication to obtain weak-interaction solutions. 
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APIENLIX A 


THE TURBULENT VORTICITY TRANSPORT EQUATION 
INCLUDING CURVATURE EFFECTS 


For flow over a curved wall a curvilinear coordinate system is introduced 
in which the x-axis is along the wall and the y-axis perpendicular to the wall. 
Hie coordinate system, therefore, consists of a set of curves parallel to the 
wall and a set of straight lines perpendicular to the wall. As shown by ref. 
15, if the curvature is denoted by k the streamwise and transverse momentum 
equations become 


it- * 7^7 & ‘ v if * l+ty * rrtyijir* 


viscous terms ( A_1 ) 
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In the curvilinear coordinate system the vorticity, Tu , is given by 


where e^ is the unit vector perpendicular to the x-y plane. The vorticity 
transport equation is obtained by subtracting Eq. (A-l) multiplied by (l-4ty) 
from Eq. (A -2) to obtain a transport equation for tiT. In the present effort 
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As in the case of Cartesian coordinates, the additional turbulence stress 
terms are obtained by dividing u and v into mean and fluctuating parts, 
averaging, and then combining the equations. Since all extra turbulence stress 
terms result from the nonlinear convective terms, the nonlinear convective 
terms are now examined in detail. Before dividing the nonlinear convective 
terms into mean and fluctuating parts, the streamwise and transverse momentum 
equations, Eqs . (A-l) and (A-2), are modified through the continuity equation 
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When the velocities are divided into mean and fluctuating parts and the mean 
continuity equation is applied the results are 
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Hie fluctuating averaged terms in Eqs. (A-8) and (A -9) contribute to the 
turbulence stresses. The contribution of the fluctuating terms to the 
vorticity equation is obtained by bringing the terms to the right-hand side of 
the equations, taking the derivative with respect to y of Eq. (A-8) multiplied 
by (l+ky), and subtracting this from the derivative with respect to x of 
Eq. (A -9) which gives the contribution 
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When the usual boundary layer type assumption is made that derivatives with 
respect to y are much large than those with respect A 5 x, this reduces to 
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By analogy to the procedure for two-dimensional flow (-u'v') is related tc the 
mean flow field through an eddy viscosity, v t . It should be noted that in the 
present calculations v-t is not solely dependent on local mean flow conditions 
but rather is dependent upon the flow hi stor y through t v :<~ turbulence kinetic 
energy equation. Hie relation between -u’v' and is given by 


- U V = -VfUt 


(A-12) 


Therefore, 
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It should be noted that the relationship between -u'v’ and given in Eq. 
(A-12) is somewhat different than the usual formulation 
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However, examining Eq. (A-3), with the assumption that Av/Ax is small it is 
found that 
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and, thus, Eq. (A-12) may be written 
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Thus it is evident that the formulation of Eq. (A-12) differs from the usual 
formulation of Eq. (A-l4) only in regions of significant curvature, or only 
in the region of the leading edge of the airfoil where the flow is expected 
to remain laminar and vt is negligible. In other regions of the airfoil, 
where the curvature is negligible, the two formulations are equivalent. 
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THEORY OF GIESING (REF 10) 



Figure 2. — Comparison between theoretically predicted and experimentally measured aerodynamic coefficients. 
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Figure 3. - Companion batwaan numarieal tolution ami analytical solution for Hiemenr flow. 
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Figure 4. - Comparison between numerical solution and analytical solution for flow about circular cylinder. 
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Figure 9. - Comparison between theoretical prediction and experimental data for 
oscillating turbulent boundary layer. 
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Figure 12. - Comparison between theoretically predicted and experimentally measured 
differential pressure coefficients for Case I . 
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Figure 13. - Comparison between theoretically predicted and experimentally measured 
differential pressure coefficients for Case I . 
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Figure 14. - Comparison between theoretically predicted and experimentally measured 
differential pressure coefficients for Case I . 
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DISTANCE FROM LEADING EDGE. X/C 


Variation of tha pranura coafficiant on tha airfoil turfaca at various initantanaoi 
angle* of attack for Casa I. 
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Figure 17 . - Variation of ikin friction coefficient in the stagnation region at various instantaneous 
angles of attack for Cate I . 
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Figure 20. - Variation of dimensionless mixing langth in tha stagnation region at various 
instantaneous angles of attack for Casa i . 
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( Figure 28 . - Variation of skin friction coefficient alo g airfoil surface in the pressure side fully turbulent 

region at various instantaneous angles of attack for Case I . 
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MOMENTUM THICKNESS, (0/0x12 



Figure 29 . - Variation of momentum thickness along airfoil surface in the pressure side fully turbulent 
region at various instantaneous angles of attack for Case I . 
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Figure 30. - Variation of displacement thick nau along airfoil Hirfaca ir> tha praMura side fully turbulent 
region at various initantanaoui angles of attack for Ca*a I . 
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Figure 31 . - Variation of dimansionlass mixing langth along airfoil turfaca in tha pressure sida fully 
turbulant region at various instantanaous anglas of attack for Casa I . 
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DISTANCE ALONG AIRFOIL SURFACE. S/C 

Figure 34. -Variation of skin friction coefficient in the s epar a ted region of the airfoil suction side at various 
in stant aneous angles of attack for Case I . 
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Figure 35.- Variation of momantum thick nan in tha saparatad tagion of tha airfoil suction sida at various 
in sta ntanaous a n alas of attack for Casa I . 
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Figure 37. - Variation of dimensionless mixing length in the separated region of the airfoil suction side at 
various instantaneous angles of attack for Case I . 
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Figure 39.- Variation of skin friction coefficient along the airfoil surface in the suction side fully turbulent 
region at various instantaneous angles of attack for Case I . 
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Variation of momentum thick nan along ttw airfoil surface in tha suction side fully turbulent 
region at various instantaneous angles of attack for Casa I . 
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DISTANCE ALONG AIRFOIL SURFACE, S/C 


Figure 42 . - Variation of dimansionlass mixing length along the airfoil surface in tha suction sida fully 
turbulent region at var i out inttantanaous angles of attack for Casa I . 
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Figure 43. - Variation of skin friction coefficient along the airfoil surface in the suction side trailing 
edge fully turbulent region at various instantaneous angles of attack for Case I . 
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Figure 45. - Variation of displacement thickness along the airfoil surface in the suction side trailing 
edge fully turbulent region at various instantaneous angles of attack for Case I . 
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Figure 47. - Predicted instantaneous velocity profiles at a fixed streamwise location 
in the suction side fully turbulent region for Case I . 
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Figure 49 . - Companion betwe en theoretically predi c ted and experimentally 
m anured differential pressure coefficients for Case II. 
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Figure 50 . — Comparison botwo a n theoretically predicted and experimentally 
meewired differential prenure coefficients for Case II . 
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ARROW INDICATES DIRECTION OF INCREASING TIME 



Figure 53. - Location of stagnation point as a function of angle of attack 
for Case II . 
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Figure 57. 


Variation of dimensionless mixing length in the stagnation region at 
various instantaneous angles of attack for Cate II . 
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Figure 58. - Variation of tkin friction coefficient along the airfoil surface in the preature tide traniition 
region at various inatantanaou* angiea of attack for Caaa II . 
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Figure 59. - Variation of momentum thickness along the airfoil surface in the pressure side transition 
region at various instantaneous angles of attack for Casa II . 
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Figure 60. - Variation of diiplaeamant thickness along the airfoil surface in the pressure side transition 
region at various instantaneous angles of attack for Case II . 
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Figure 63. - Variation of momentum thick ness along tha airfoil surface in the pressure side 
fully turbulent region at various instantaneous angles of attack for Case II . 
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Figure 65. 


Variation of dimensionless mixing length along the airfoil surface in the pressure 
side fully turbulent region at various instantaneous angles of attack for Case II . 
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Figure 69. - Variation of displacement thickness in the separated region of the airfoil suction side at various 
instantaneous angles of attack for Case II . 
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Figure 70. - Variation of dimensionless mixing length in the seperated region of the airfoil suction side at 
vario instantaneous angles of attack for Case II . 
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Figure 72. - Variation of momentum thick neia along tha airfoil surface in tha auction sida 
iully turbulent ragion at variout instantanaoua angles of attack for Caw II . 





Figure 73. - Variation of displacement thickness along the airfoil surface in the suction side 
fully turbulent region at various instantaneous angles of attack for Case 1 1 . 
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Figure 74. - Variation of dimensionless mixing langth along tha airfoil surfsc? in tha auction side 
fully turbulent region at various inatantanooui angles of attack for Casa i I . 
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Figura 76. - Variation of momentum thickness along the airfoil surface 
in the suction side trailing edge fully turbulent region at 
various In s tantaneous angles of attack for Casa 1 1. 
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Figure 77. Variation of displacement thickness along the airfoil surface in the suction side 
trailing edge fully turbulent region at various instantaneous angles of attack 
for Case 1 1 . 
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THEORY OF GIESING (REF. 101 


EXPERIMENTAL DATA OF GRAY AND LIIVA IREF 31 
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Figure 81. - Comparison between theoretically predicted and 

experimentally me ssured aerodynamic coefficients 
for Case III . 
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Figure 84. — Location of statnation point as a function of normalized transverse location 
for Case III . 
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Figura 88. - Variation of momantum thick nata in tht ttaonation region at 

variout inatantanaouf nonnaliiadtra na vana locations fot Caw III . 
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Figure 87. - Variation of ditplacemant thickness in the stagnation region at 
various instantaneous normalized transverse locations for Case ill . 
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Figure 88. - Variation of dimanslonlass mixing langth in tha stagnation ragion at various 
instantanaous normalizad transvarsa locations for Casa III . 
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Variation of momentum thickness along the airfoil surface in the pressure side 
transition region at various instantaneous normalized transverse locations 
for Case III . 
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Figure 91. - Variation of displacement thickness along the airfoil surface 
in the pressure side transition region at various instantaneous 
normalized transverse locations for Case III . 
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Figure 92. - Variation of dimension lass mixing length along the airfoil surface in 
the pressure side transition region at various i ns tantaneous normalized 
tr a nsverse locations for Casa III . 
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Figure 93. — Variation of skin friction coefficient along the airfoil surface 
in the pressure side fully turbulent region at various 
instantaneous normalized transverse locations for Case III . 
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Figure 94. — Variation of momentum thickness along the airfoil surface in the pressure tide 
fully turbulent region et verious instantaneous normalized transverse locations 
for Case ill. 
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Figure 100. - Variation of ditpleoamant thick nan in tha reparation region of tha airfoil 
auction aida at varioui inatantanaoua normaiizad trantvarre Iccationi 
for Care III . 
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Figure 103. - Variation of momentum thickness along the airfoil surface in the suction side 
fully turbulent region at various instantaneous normalized transverse locations 
for Case III . 
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